8.27 Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

8.27 Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

8.27 Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra L Benda, P Sˇteˇpa´nek, J Kaminsky´, and P Bourˇ, Institute of Organic Chemi...

1MB Sizes 78 Downloads 269 Views

8.27 Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra L Benda, P Sˇteˇpa´nek, J Kaminsky´, and P Bourˇ, Institute of Organic Chemistry and Biochemistry, Prague, Czech Republic r 2012 Elsevier Ltd. All rights reserved.

8.27.1 8.27.2 8.27.2.1 8.27.2.2 8.27.3 8.27.3.1 8.27.3.2 8.27.4 8.27.4.1 8.27.4.2 8.27.4.3 8.27.4.4 8.27.4.5 8.27.5 8.27.5.1 8.27.5.2 8.27.6 8.27.6.1 8.27.6.2 8.27.6.3 8.27.7 8.27.7.1 8.27.7.2 8.27.7.3 8.27.8 8.27.9 References

Introduction: Chiral Techniques and Interpretations of the Spectra Overview of Ab Initio and DFT Methods Theory of Quantum-Chemical Methods Quantum-Chemical Calculations of Spectroscopic Properties Optical Rotatory Dispersion Definition Applied Computations Electronic Circular Dichroism Experimental Quantities Simplified Models First-Principles ECD Computations Origin Dependence of Theoretical ECD Applications Magnetic Circular Dichroism as a Chiral Method for Achiral Molecules Theory of MCD Practical Calculations Vibrational Circular Dichroism Experimental Aspects VCD Theory Spectral Simulations Raman Optical Activity History and Experiment ROA Theory Simulation Examples Common Aspects of the Computational Approaches Conclusions

Glossary Chiroptical spectroscopy It explores differences in absorption, scattering, or dispersion of left- and rightcircularly polarized light on molecules. Density functional theory (DFT) It is an efficient method for computation of molecular energies, based on electronic density. Electronic circular dichroism (ECD, also CD or UVCD) It is the difference in absorption of left- and rightcircularly polarized light by a chiral molecule (which has neither plane nor center of symmetry). Magnetic circular dichroism (MCD) It is circular dichroism in the presence of a static magnetic field.

8.27.1

520 522 522 524 525 525 526 527 527 528 528 529 530 531 532 533 533 533 533 535 536 536 536 538 539 542 542

Optical rotatory dispersion (ORD) It is a manifestation of the different speed of left- and right-circularly polarized light, measured also as a rotation of linearly polarized light. Raman optical activity (ROA) It is a difference in Raman scattering of right- and left-circularly polarized light. Response theory It is a general formulation of the timedependent perturbation theory (including DFT) allowing to calculate molecular properties. Vibrational circular dichroism (VCD) It is circular dichroism in the infrared spectrum, based on vibrational transitions in molecules. Vibrational optical activity (VOA) It usually involves VCD and ROA.

Introduction: Chiral Techniques and Interpretations of the Spectra

According to a classical wisdom, molecules are too small to be seen by the naked eye because their size is smaller than the light wavelength. This is certainly true if we think of them as kinds of sand grains or small beetles. However, mere seeing may not be the most important point for us; the knowledge acquired from it is. In this respect, the many spectroscopic techniques provide us with much more information about molecular behavior, such as their flexibility, structure, and interactions, than it would be possible by visualizing molecular shape, although, to be fair, we must admit that other techniques, such as tunneling microscopy, also allow us to observe the shape nowadays.

520

Comprehensive Chirality, Volume 8

http://dx.doi.org/10.1016/B978-0-08-095167-6.00848-X

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

521

It is also true that NMR (see Chapter 8.28) and X-ray spectroscopy (see Chapter 9.16) provide higher spatial resolution than the optical spectra. For practical applications, it is always important to realize the advantages and drawbacks of the plethora of methods available. The low-resolution experiments are usually cheaper, faster to perform, and various parameters (pH, temperature, concentration) can be changed more easily. For many compounds that do not crystallize, X-ray and similar diffraction techniques cannot be used. A very typical application area for optical spectroscopy is the monitoring of conformational changes; for a large protein, for example, the all-atom or high-resolution information would be redundant if we only look at a variance of the secondary structure. In some cases, even conformers not resolvable by NMR can be discriminated by the optical methods. This is actually quite natural, as the optical response (absorption/emission of a photon) is very fast and the spectrum measured is an algebraic sum of the subspectra. However, the NMR technique is much slower and often records only an average signal from many molecular forms. The chirality aspect of optical polarized spectroscopies is an important advantage. The conventional NMR is blind to chirality.1 Needless to say, following the chirality along synthetic pathways is important for organic chemists, as it is for quality control in the drug and food industry (see Chapters 1.1 and 1.7). Monitoring the different interactions of the left- and right-circularly polarized light is also generally useful, as the spectra become more sensitive to molecular shape and finer conformational changes than for unpolarized techniques. As we will see, this is caused by the nature of the chiral phenomena: whereas the nonpolarized spectral response (e.g., absorption) is governed by the electric dipolar interaction, both electric and magnetic dipole moments are needed for the chiral signal (e.g., circular dichroism). To be able to utilize the chiral spectra, we have to understand them. From the beginning, the experiment was coupled with theory. An observation of the optical rotation (OR) had to be related to the molecular level by Louis Pasteur in 18472 to be useful for chemists. Also, later, the chiral optical phenomena intertwined with physics, including the theory of quantum electrodynamics3 and the parity violation observation.4 Not all fundamental problems have been resolved as yet; for example, according to the Schro¨dinger equation, the left- and right-hand enantiomers are strictly speaking only temporary forms that must spontaneously interconvert, which has, however, never been observed for most compounds. The parity violation causes heavy atoms to be optically active;5 these aspects, however, are beyond the scope of our overview. For chemical applications, we stay with the classical concept of chiral ‘mirror symmetry’ molecules or samples, probed by the circularly polarized light (Figure 1). The phenomenon can easily be understood by realizing that the left- and right-hand (or screw) symmetry of the sample causes different interactions with the light. Because of the duality principle, the light can be thought of either as a photon spinning in one or the other direction along the course of propagation or as an electromagnetic wave with the field vector (perpendicular to the propagation) rotating as a left or a right screw. For ‘achiral’ molecules, that is, those that have a center or a plane of symmetry, chiral phenomena can be induced by an external magnetic field.

L

R

(l R − l L) × 103

20

L-alanyl-L-alanine

10 0 −10 −20

D-alanyl-D-alanine

80

lR + lL

60 40 20

L-alanyl-L-alanine

D-alanyl-D-alanine

0 200

400

600

800

1000 1200 1400 1600 1800

Wavenumber (cm−1) Figure 1 Spectroscopic exploration of molecular chirality: (þ) and () isomers of Ala–Ala react differently with the circularly polarized light and provide opposite Raman Optical Activity spectra (red and blue). This cannot be achieved through an unpolarized measurement, where the isomers provide the same Raman spectrum (black). The spectra were recorded on a spectrometer located at the Charles University, Prague.

522

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Also, we exploit the so-called semiclassical approximation. This means that the molecule is treated by the rigorous method of quantum mechanics, whereas the incident light is treated classically. As an electromagnetic wave, the radiation inflicts perturbations upon the quantum system. In other words, we do not care about the quantum properties of the light. Even with these simplifications, there were many obstacles that had to be overcome before the spectral simulations could be used routinely. An excellent example is the electronic contribution to the vibrational circular dichroism that disappears in the Born–Oppenheimer approximation. This led to emphatic discussions among the scientists, as they did observe the spectra readily in their laboratories. Finally, using algebra, Stephens showed that this phenomenon can also be determined using the general principles of quantum computational chemistry.6 Another common problem of chiral methods that needs to be mentioned is the well-known dependence of magnetic properties determining the chiral molecular response on the choice of origin if calculated using approximate methods. Today, the use of the gauge including (or independent or invariant) atomic orbitals (GIAO) is considered to be the most universal approach to eliminate origin dependence. Because of the dependence only on the electronic ground state, but also because of the better resolution in the vibrational region, molecular vibrational optical activity (VOA, namely Raman optical activity, ROA, and vibrational circular dichroism, VCD) can usually be simulated more reliably than optical activity involving electronic transitions (electronic circular dichroism, ECD, magnetic circular dichroism, MCD, optical rotatory dispersion, ORD). In any case, available computer programs provide many means how to support all these spectroscopic techniques. The possibility of consistently simulating the spectra has led to the commercial success of instrument manufactures: whereas the first attempts to commercialize VCD in the 1980s failed, after 1997, the technique successfully reappeared on the market. Soon, the first ROA spectrometers were sold as well. Nevertheless, the development of individual techniques may occur in the future. For example, because it has not been possible to simulate MCD reliably until B2008, it may be interesting to see whether the availability in common computer codes also stimulates further development of this rather rare technique. Still, more accurate simulations, taking into account molecular flexibility and interactions with the solvent, are needed to make the chiral spectroscopies more useful in molecular structural studies. Before discussing detailed computational aspects, we should mention that technological advances are also expected and needed (e.g., the time-dependent VOA based on novel detection schemes,7 higher sensitivity of VOA instruments, resonance and surface enhancement of the Raman techniques). The discoveries of new phenomena, such as the optical activity enhancement in macromolecular-ordered structures observed in the vibrational region,8 cannot be excluded. The interplay between experiment and theory thus governs the development of the techniques. All these developments indicate that chiral spectroscopy techniques combined with theoretical modeling will remain a dynamical and useful means for studying molecular structures and properties.

8.27.2

Overview of Ab Initio and DFT Methods

The methods of quantum chemistry represent the most rigorous but also very demanding approach for theoretical predictions of chiroptical spectra. In practice, we can distinguish between so-called wavefunction (or ab initio) methods and methods based on the electron density. Note, however, that in theory, both approaches are equivalent. Calculations can yield band positions (vibrational frequencies for VCD and ROA, electronic excitation wavelengths for ECD, MCD, and ORD) as well as spectral intensities and many other underlying molecular properties. For vibrational spectra (VCD and ROA), the ground-state molecular properties are conveniently expressed as mixed derivatives of the total electronic energy E and electromagnetic properties with respect to the electric field E, electric field gradient rE, magnetic field B, nuclear displacements Rl, and momenta Pl; for the electronic properties, excited state wave-functions and transition integrals are needed (Table 1).

8.27.2.1

Theory of Quantum-Chemical Methods

There are many methods more or less approximating the Schro¨dinger equation HˆC ¼ EC that can be used for evaluating the total electronic energy E and its derivatives. The starting point for wavefunction ab initio methods is the Hartree–Fock (HF) approximation, where the Ansatz of the Slater determinant is imposed on the multidimensional wavefunction C. Simpler one-particle HF HF spin-orbitals cHF a are found together with orbital energies ea as solutions of the HF equations for N electrons " # N  2 X  ^f HF cHF ¼  _ r2 þ vðrÞ þ ^J b  K ^ b cHF ¼ eHF cHF ; ð1Þ a a a a 2me b¼1 HF where _ is the reduced Planck constant and me is the electron mass. The Fock operator ^f consists of the electron kinetic energy ^ b . HF calculations are term, the one-electron nuclear potential v(r), and two-electron Coulomb and exchange operators ^Jb and K complicated mainly due to the presence of two-electron operators that implicitly depend on the HF solutions cHF a . The HF equations are thus nonlocal and must be solved iteratively until convergence. Each spin-orbital cHF a is a product of spatial and spin parts c(r,s) ¼ j(r)s(s), where j(r) is the molecular orbital (MO) and s(s) is a spin function. In the most common closed-shell molecular systems, two electrons of opposite spin share the same MO.

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Table 1

523

Molecular properties used for OA simulations

Property Vibrational domain Force field (FF), Hessian, second energy derivatives Dipole derivatives, atomic polar tensor (APT) Magnetic dipole derivatives, atomic axial tensor (AAT) Polarizability derivatives G0 -polarizability derivatives A-polarizability derivatives Electronic spectra Dipole strength

Definition

Needed for

q 2E q Ral q Rbn q ma q Rel q ma q Pel q aab q Rel q G 0 ab q Rel q Aa;bg q Rel

Vibrational frequencies

 hg j^ ljk ij2 ^ jg i Imhg j^ ljk ihk jm Equation 28 Equations 12 and 42

Rotational strength A, B, C terms G0 -tensor

MOs are typically obtained as linear combinations of atomic orbitals fm (AO), X cma jm ; ja ¼

IR intensities, VCD VCD Raman and ROA intensities ROA intensities ROA intensities

Absorption intensities ECD MCD ORD

ð2Þ

m

where cma are expansion coefficients. The usual form of AOs is based on electron orbitals derived for hydrogen-like atoms, which, for this purpose, are simplified and combined into real functions. Let us recall that the exact hydrogen orbitals can be written as fnlm(r, W, j) ¼ Rnl(r)Ylm(W, j), where n, l, and m are orbital quantum numbers, the angular parts Ylm(W, j) are spherical harmonics, and the radial parts Rnl(r) ¼ Pnl(r)ear are products of Laguerre polynomials Pnl(r) and the exponential function. Basis sets in quantum-chemistry calculations differ the most in the AO radial part. Basis sets of the form Rnl(r)Bear are called Slater-type orbitals (STO). However, calculation of the two-electron integrals with STOs is relatively complicated; of the major quantum-chemistry packages, only ADF uses STOs. Other implementations circumvent the use of STOs by fitting the AO radial P 2 part with linear combinations of Gaussian functions, Rnl ðrÞEPnl ðrÞ n cn ean r , which yields Gaussian-type orbitals (GTO). The advantage of GTOs is the possibility of rapid evaluation of the two-electron integrals, which mostly compensates for the need for a larger number of Gaussians compared with STOs. The HF ground-state energy EHF ¼

N X a¼1

eHF a 

N X N     1X ^ b cHF cHF ^J b  K a 2 a¼1 b¼1 a

ð3Þ

already provides a good estimate of the exact Schro¨dinger ground-state energy E. The remaining error, the correlation energy Ecorr ¼ E  EHF, may be approximated using variational or perturbational methods. Most often, the Møller–Plesset (MP) perturbation theory and the coupled-clusters (CC) techniques are used. However, correlated ab initio calculations of chiroptical properties are still relatively rare, and are not covered much in this brief overview. An alternative solution to the many-electron problem is provided by the density functional theory (DFT).9 It focuses on the electron density r(r) that is (according to the Hohenberg–Kohn theorem) uniquely related to the one-electron potential v(r) and thus unambiguously determined by nuclear positions, similar to the wavefunction. The electronic energy is then defined as a functional (function of a function) of the electron density, E ¼ E[r(r)]. In a very convenient DFT formulation by Kohn and Sham,10 the stationary condition for energy minimum dE[r(r)]/dr(r) ¼ 0 is equivalent to a set of self-consistent Kohn–Sham (KS) equations   Z 2 rðrÞ KS KS ^f KS cKS ¼  _ r2 þ vðrÞ þ 1 dr0 þ vxc ðrÞ cKS a a ¼ ea ca 0 4pe0 jr  r j 2me rðrÞ ¼

N  X  cKS ðrÞ2 a a¼1

ð4aÞ

ð4bÞ

524

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra vxc ðrÞ ¼ dExc ½rðrÞ=drðrÞ

ð4cÞ

for KS orbitals cKS a ; vxc(r) is the exchange-correlation potential. The orbitals are, strictly speaking, not needed in DFT, but without them it would be very difficult to calculate the kinetic part of the energy precisely. Although the exact form of the exchangecorrelation functional Exc[r(r)] is not known, a large number of approximate DFT functionals are available today. It should be emphasized that the performance of different DFT functionals varies and should be tested (Table 2). The KS ground-state energy is

EKS ¼

N X a¼1

Table 2

eKS a 

1 1 2 4pe0

ZZ

rðrÞrðr0 Þ 0 dr dr  jr  r 0 j

Z

vxc ðrÞrðrÞdr þ Exc ½r:

Overview of DFT functionals

Basic types Pure local Pure nonlocal, GGA Hybrid DFT-D Double-Hybrid

Dependent on density (r) Dependent on density derivatives (rr,y) Mixed with HF exchange Explicitly corrected for dispersion (van der Waals) energy Mixed with MP2 correlation

Our experience Good for force field and vibrational spectra Good for CD Good for MCD

B3LYP, B3PW91, BPW91 (n(C¼O) region) B3LYP, CAM-B3LYP SAOP, B3LYP

Not recommended for FF Not recommended for CD

8.27.2.2

ð5Þ

HF, B1LYP, PBELYP, PW91LYP, BP86, BLYP, G96LYP, PW91PW91, SVWN5, LSDA, OLYP, VSXC, HCTH, M5, M06 BP86, BPW91, BLYP (GGA in general)

Quantum-Chemical Calculations of Spectroscopic Properties

Derivatives of the ground-state energy needed for the spectroscopic properties can be obtained in several ways. The simplest is the numerical differentiation, where the energy derivatives are calculated from finite energy differences, for example,

q 2E EðDx; DyÞ  EðDx; DyÞ  EðDx; DyÞ þ EðDx; DyÞ ; E 4DxDy q xq y

ð6Þ

where Dx and Dy are nuclear displacements (typically B0.02 A˚). The numerical schemes are simple but have limited accuracy, and often lead to tediously long computations. Numerical differentiations can also be used in combination with analytical techniques. For example, we can numerically obtain second energy derivatives from analytical first derivatives (gradients). Analytical differentiation techniques are more convenient for reliable calculations. Using the time-independent perturbation theory, for example, the energy derivatives can be expanded as sums over excited states 9CkS ¼ 9kS. The second energy derivatives for a system with a nondegenerate ground state 9CgS ¼ 9gS can then be expressed as X hgj qqH^x jkihkj qqH^y jgi ^ q 2E q 2H ¼ hg j jgi þ 2Re : q xq y q xq y Eg  Ek ka g

ð7Þ

Such sum-over-states (SOS) expansions are convenient for gaining a physical insight into spectral parameters. However, the calculations with approximate excited states can be rather inaccurate. An alternative approach is provided by the coupled-perturbed (CP) techniques.11 The HF or KS equations can be rewritten for a molecular system in the presence of an external perturbation, such as an electric or a magnetic field. Thus, we obtain CP equations ð1Þ ð0Þ ð^f  eða1Þ Þjða0Þ  ð^f  eða0Þ Þjða1Þ ¼ 0

ð8Þ

ð1Þ ð1Þ ð1Þ comprising the first-order orbital and energy corrections ja and ea , the first-order Fock operator ^f , and the unperturbed ð0Þ ð0Þ ð0Þ ð1Þ quantities ^f , ja , and ea . In the usual formulation, the first-order orbital corrections ja are expanded into unperturbed

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

525

orbitals jða1Þ ¼

  E D ð0Þ  ð1Þ  ð0Þ X jr ^f ja ð 0Þ

ð0Þ

ea  er

r

jðr0Þ

ð1Þ

ð9Þ ð1Þ

In both the HF and the KS cases, the first-order operator ^f depends on the perturbed orbitals ja . The CP equations are thus ð1Þ solved iteratively. According to the ‘2n þ 1’ rule, the first-order orbital corrections ja are sufficient to evaluate the second-order energy derivatives. Analytical expressions for higher order energy derivatives are obtained analogously. The time-independent perturbation theory provides molecular properties at the limit of static external perturbations. Using this approach, harmonic frequencies, IR/VCD intensities, and nonresonance Raman/ROA intensities can be obtained. Other spectroscopic properties such as excitation energies, ECD, MCD, and ORD intensities and frequency-dependent Raman/ROA intensities have to be treated via the time-dependent perturbation theory. One of its most general formulations is the response theory.12 Most of the spectroscopic properties are correctly described already at the lowest (linear) level of response theory. Special cases of the linear response theory appeared earlier as the time-dependent HF theory13 (TDHF; also referred to as the random phase approximation, RPA) and the time-dependent DFT14 (TDDFT). The response theory provides excitation energies and spectroscopic properties when the linear response equations are solved. However, the excited states are not needed explicitly. For example, calculation of the dynamic polarizability by perturbation expression,

X hgj^ ma jkihkj^ mb jgi hgj^ ma jkihkj^ mb jgi þ ; ð10Þ aab ðoÞ ¼ Ek  Eg  _ o Ek  Eg þ _ o ka g ^a are components of the requires exact excited states 9kS; m dipolemoment operator. In the response theory, the polarizability can ^a ; m ^b o that can be accurately evaluated using perturbed ground-state be identified with a linear response function, aab ðoÞ ¼  m HF or KS molecular orbital derivatives. Theoretical modeling of chiroptical spectra requires the calculation of magnetic molecular properties (such as AAT or the G0 -tensor). In a complete basis-set limit, observable magnetic properties are independent of the choice of common origin. However, finite basis sets can introduce significant origin-dependent errors that need to be removed for reliable calculations. The most frequent way to remedy the ‘gauge’ problem is the introduction of the gauge-including atomic orbitals, GIAOs15 (or also London orbitals) h e i ðrÞ ¼ exp i Al  r fl ðrÞ ð11Þ fGIAO l _ where the vector potential Al ¼ B  Rl/2 at nucleus l is explicitly included in the complex prefactor. The response equations must then be solved with the GIAOs, that is, with MOs dependent on the magnetic field. Technically, however, the same AOs (basis functions) are input in actual GIAO calculations as for the origin-untreated case. It is often claimed that magnetic molecular properties calculated using GIAOs (ROA, VCD, ECD (in some formulations), but also magnetic susceptibilities and NMR shielding) converge faster with the size of atomic basis set than for other procedures eliminating the origin dependence.

8.27.3

Optical Rotatory Dispersion

OR is a rotation of the plane of linearly polarized light when passing through the sample. It is the oldest chiral phenomenon, observed in 1815. We explain it by different speeds of left- and right-circularly polarized light components (of the linearly polarized light) in a chiral environment. For the wavelength dependence of OR, ORD is of interest. Conventionally, OR was associated with configurations of molecular chiral centers. Only the latest computational analyses have revealed that measured values often include contributions from many conformers, which can even have different signs.16,17 Thus, precise knowledge of conformer populations is necessary to reproduce the experiment,18 and vice versa, calculated OR values can be used to explore the conformational space. Since OR can be measured easily, its accurate simulations can be handy tools for structural studies (see Chapter 8.21).

8.27.3.1

Definition

The rotation is caused by molecular polarization, and both the electric and the magnetic field components of the electromagnetic radiation contribute to this phenomenon. The quantum theory of OR was initiated by the work of Rosenfeld in 1928.19 OR calculations using first-principle methods became more feasible since the work of Amos,20 who used the HF theory in a static limit as implemented in the CADPAC program.21 Many applications soon followed.22 The observed rotation is usually reported as the specific rotation [a]l in degree (dm g cm3)1, that is, normalized to unit density (g cm3) and cell length (dm) at the wavelength l (corresponding to the angular frequency o). It can be expressed as ½al ¼

7200NA o2 bðoÞ; c2 M

ð12Þ

526

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

where bðoÞ ¼ ð1=3oÞ G0aa, M is the molar mass (g mol1), NA is Avogadro’s number, c is the speed of light (m s1), and o is in rad s1. The OR parameter b(o) (cm4) is thus related to the trace of the electric dipole-magnetic dipole polarizability G0 ab defined below (equation 42). For b(o) in atomic units and o in cm1 we obtain [a]l ¼ 1.3423 104o2b(o)/M.

8.27.3.2

Applied Computations

For flexible molecules, it is desirable to find the most important minima on the potential energy surface. Depending on the complexity of the system, we can use ab initio force fields, semiempirical methods, automatic conformational searches, etc. Sometimes, the explicit involvement of solvent molecules is required, using molecular dynamics (MD) or Monte Carlo (MC) sampling. Populations of individual conformers can be obtained from the Boltzmann distribution. Apart from the conformational problem, one has to carefully consider the theoretical level (basis set, electron correlation treatment/density functionals), preferably including environmental (solvent) and vibrational corrections. OR dependence on the basis set was studied at various approximations.23–25 It is important to use a balanced set of polarization and diffuse functions. A simple addition of the diffuse functions to a small basis sets may not improve the results. We found that the aug-cc-pVDZ basis set was an acceptable compromise between accuracy and computational cost.26 The Sadlej–pVTZ basis set27 for methyloxirane yielded similar OR as a larger aug-cc-pVTZ. For small compounds, the aug-cc-pVQZ results were close to the complete basis set limit (CBS);23,28 the limit could be conveniently estimated using a simple two-point scheme.29 Chiroptical properties are also very sensitive to approximations performed in the electronic structure calculations. The basic HF theory is generally not recommended for OR. TDDFT appears to be a better alternative. It uses the origin-independent formulation,30 the computations are reasonably fast, and a reasonable accuracy can be achieved, allowing meaningful comparisons with experiment.25,31 Only sometimes TDDFT produces worse results than HF (see Figure 2),26 and rarely (e.g., anionic proline) it fails completely.17 Analysis of several aminoacids in different ionic states indicated that the errors may be caused by underestimation of electronic excitation energies.32 Hybrid functionals with a higher fraction of exact exchange (e.g., BHandH) are often more accurate than the GGA ones. Also, the CAM-B3LYP functional with a long-range correction for asymptotic behavior seems promising for OR.33

H O

H

OMe

Empirical

200

O

B3LYP0.8HF

OPBE

OLYP

B3PW91

B3LYP

300

HF

Deg (dm g cm−3)−1

H

100

0 Theoretical level Figure 2 Optical rotation [a]589 of a sugar derivative as calculated with different theoretical levels and the aug-cc-pVDZ basis set. The B3LYP0.8HF column represents modified B3LYP functional with 80% HF exchange. Empirical (experimental) values were obtained from the raw experiment by subtracting the vibrational corrections.

The high computational costs of wavefunction electron correlation methods (CI, MRSCF, CC) make them still impractical for larger and flexible molecules, even though their accuracy is higher than that for DFT and can be controlled beforehand.28,31 So far, these methods have mostly been used for benchmark calculations. The use of a solvent significantly influences OR. We can approximately attribute a solvent’s effect to two factors: influence on the conformation and change of the OR magnitude for a fixed geometry. The standard approach treats the solvent as a polarizable continuum,34 which usually improves the agreement with experiment compared with the vacuum calculations. But the polarizable continuum models (PCM) only partly explain a solvent’s influence, especially polar solvents. In particular, the effects of hydrogen bonds are difficult to describe by PCM. A simple point-charge solvent (water) representation is a computationally very efficient alternative.35 An explicit cluster averaging based on an MD simulation36 is in general more appropriate, but obviously also more computationally demanding. A strong effect of aggregation was observed for a dimer formation of pantolactone.37 Theoretical results were compared with experimental OR dependent on the pantolactone/CCl4 ratio. Beratan et al. even introduced a concept of chiral imprinting effects on the surrounding medium, which can even dominate the chiroptical response of a solute.38

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

527

In general, averaging of OR values over higher frequency vibrational molecular motions is also desirable. It involves zero-point (dominant) and finite-temperature (usually small) contributions. However, because of the limited accuracy of the computations, the addition of the vibrational corrections has not convincingly improved computed OR data so far, even though the corrections were a substantial part of the total values.26,39,40 Considering the influence of the solvent on the zero-point vibrational corrections to OR for S-methyloxirane led to limited success.40 The effect of the solvent on the vibrational corrections was approximately constant within the entire frequency range, whereas the electronic equilibrium component decreased with increasing frequency. In some solvents, the vibrational correction exceeded the equilibrium electronic part. In this respect, cavity ring-down polarimetry, enabling measurements of OR in vacuum, seems to be a powerful tool for solute–solvent interaction studies, as it provides benchmark vacuum values (see Chapter 8.21).41

8.27.4 8.27.4.1

Electronic Circular Dichroism Experimental Quantities

ECD (also called CD or UVCD) is a very popular technique for structural analysis (see Chapters 8.22, 8.23, and 8.24). The experiment is simple and does not require high concentrations. It provides richer information about the electronic transitions than standard UV–vis absorption spectroscopy: unlike that for absorption, ECD bands can be both positive and negative. However, it is applicable only to chiral molecules. The differential absorption between left- and right-circularly polarized light is DA ¼ AL  AR :

ð13Þ

A molecule is characterized by a normalized quantity, the differential molar absorption coefficient De ¼ eL  eR ¼ DA=cl;

ð14Þ

where c is the molar concentration (mol L1) and l is the cell length (by default given in cm for CD). For historical reasons, CD is also reported as ellipticity Y (in degrees) and molar ellipticity [Y] (in deg L mol1 cm1). The conversion to De is ½Y ¼ 3298 De:

ð15Þ

The integral (area) of an absorption band corresponding to the g - k electronic transition is proportional to the dipole strength D, defined as D ¼ hgj^ ljkihkj^ ljgi:

ð16Þ

^ jg i R ¼ Imhgj^ ljkihkjm

ð17Þ

Similarly, in ECD, the rotational strength

is introduced so that D ¼ kD

Z

eðlÞ dl; l

ð18Þ

DeðlÞ dl: l

ð19Þ

ABSband

and R ¼ kR

Z CDband

In SI units, kD ¼ 1.022 1060 and kR ¼ 7.659 1054. When the strengths are in debyes, kD ¼ 9.184 103 and kD ¼ 2.296 103. P According to a sum rule42 ( i Ri ¼ 0, summing over all transitions), the CD spectrum is conservative within a large interval of wavelengths. The rotational strength is a pseudoscalar, and it changes its sign upon coordinate inversion. Sometimes, a dissymmetry factor is introduced as g ¼ De=e ¼ DA=A:

ð20Þ

For usual systems, g is of the order of 103. We can see that a molecule has nonzero CD if the scalar product of the transition electric and magnetic transition dipole moments (17) does not disappear. This means that the molecule must absorb the light, the magnetic moment must not be zero, and the electric and magnetic moments must not be orthogonal.

528

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

ECD and ORD spectra can be mutually converted via Kramers–Kronig transformations.43 For example,

½aðlÞ ¼

2 p

ZN 0

½Yðl0 Þ

l0 dl0 l2  l02

ð21Þ

if both molar rotation [a] and molar ellipticity [Y] are in the same units. For this purpose, several numerical methods have been developed.44

8.27.4.2

Simplified Models

Some molecules can be approximately thought of as ensembles of individual chromophores. For example, the nucleic acid bases or the amide groups in peptides are only weakly interacting with each other and with other molecular parts. Electronic transitions can thus be calculated separately for individual chromophores at a quantum-mechanical level. Even more often, transition dipoles and energies are considered as empirical parameters. The interaction between the chromophores can be treated in a simplified way. The coupled oscillator (CO) method, also referred to as the transition dipole coupling (TDC), approximates the interactions by the dipole–dipole term. One way of deriving the TDC expressions is by using deVoe’s polarizability theory.45 We provide a simplified derivation: let a ground state of a system of M chromophores be a product, 9E0S ¼ g1g2? gi? gM. Singly excited states 9EiS are obtained from mono-excitations on any chromophore as 9EiS ¼ g1g2? gi? gM, i ¼ 1,y, M. The individual chromophore states yield the Schro¨dinger equation Hˆigi ¼ egigi (the energy scale is chosen so that egi are zero) and Hˆiki ¼ egiki. As shown above, the energies ei and transition dipole moments li ¼ hgi j^ ljki i can either be calculated quantum mechanically for each chromophore or taken as empirical parameters. For identical chromophores, we use the degenerate CO model (DCO, all ei are the same). P ^ i þ V; ^ V ^ is the interaction potential. The total wavefunction for a state l is a ^ ¼ The Hamiltonian of the entire system is H iH P l c linear combination of the monoexcitations jCl i ¼ M j i¼1 i Ei i; the ground state does not couple with the excited states in this approximation. Then we can write the Schro¨dinger equation in the matrix form: Hij clj ¼ el c li . The diagonal elements Hii ¼ ei, and the off-diagonal elements are approximated by the dipole–dipole interaction Vij ¼ (4pe0)1 rij2 li  lj  3li  rij lj  rij =rij5 . For each P ljWl i ¼ j clj lj . For the magnetic dipole moment, we realize its transition g-k, we obtain the transition electric dipole as hE0 j^ origin dependence: each electric dipole positioned at ri contributes to the magnetic dipole by mi ¼ ð1=2Þri  q li =q t, so that P ^ jWl i ¼ i=2_ j clj ej rj  lj , where we used q li =q t ¼ iðei =_ Þli . hE0 jm The transition moments obtained are used in equations 16 and 17 to calculate the absorption and CD intensities. Of many computer implementations, we can consider, for example, the CO program developed by Hug et al.46 The CO scheme can be developed into a more robust matrix method.47 The system wavefunction is computed in the same way as for CO on the basis of individual chromophores, and the transition energies are then obtained by the diagonalization of a Hamiltonian matrix. The Hamiltonian, however, can be systematically improved, including accounting for a charge transfer. The matrix method has been used, for example, to analyze CD spectra of proteins.48 On the webpage of Hirst’s group, one can find an online program called DichroCalc49 for the calculation of CD spectra of proteins based on the matrix method (http://comp. chem.nottingham.ac.uk/dichrocalc/).

8.27.4.3

First-Principles ECD Computations

First-principles ECD simulations are currently feasible for fairly large molecules. The bottleneck is the excited state description; the rotational strengths require only a small additional effort. In general, it is necessary to take into account the electron correlation for reliable results. Large basis sets are also needed; however, according to our experience, split valence basis sets with diffuse functions (e.g., 6–31 þ þ G) already provide satisfactory results for many systems. Before focusing on ECD, let us at least briefly recall the main methods that determine electronic excited states. They can roughly be divided into three groups (Table 3): post-HF or so-called wavefunction methods, computations based on the timedependent KS or HF equations, and semiempirical approaches. In the configuration interaction (CI) approach, the ground-state HF wavefunction is enhanced by a linear combination of other Slater determinants. These are obtained from HF ground-state orbitals, with various electronic configurations in occupied and virtual orbital spaces. If all possible determinants and infinite basis sets could be included, we can consider the full CI (FCI), providing an exact solution of the Schro¨dinger equation. Obviously, a truncation is needed in practice. CI with single ‘excitations’ is called CIS. Typically, CIS excitation energies are overestimated due to the lack of a correlation energy. The CIS method is variational and size-consistent. The complete active space self-consistent field (CASSCF)57 is a much more advanced approach than CIS. However, its computational cost increases exponentially with the size of the active space. Moreover, the choice of the active orbital space is not trivial and must be adapted to a particular system. If it is too small, a part of the correlation may be neglected, which results in a substantial error. The method can be combined with perturbation computations in the CASPT252 variant.

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Table 3

529

Common methods for excited state calculations

Theoretical basis

Methods

Wavefunction methods

FCI, CIS50, CISD CASSCF51 CASPT252 EOM-CCSD53,54 SAC-CI, CC2, CC35

DFT

TDDFT, TDHF/RPA55

Semiempirical

INDO/S, ZINDO56

An efficient treatment of the configurations for single-reference problems is based on the coupled cluster (CC) algorithm.58 A widely used CCSD scheme includes single and double excitations. The CC theory was optimized for the calculation of excited states via an equation of motion or linear response theory.54 A method similar to CC is the symmetry-adapted cluster CI (SACCI)59 and other CC schemes of second and third order (CC2, CC3).60 Traditional quantum-chemical semiempirical methods consist of various simplifications in calculations of the two-electron integrals, such as neglect of differential overlap. Ionization potentials can be used as parameters. The common CNDO/S (complete neglect of differential overlap-spectroscopic parameterization) and ZINDO/S (Zerner’s intermediate neglect of differential overlap) schemes are based on the CIS approach. The semiempirical methods are very fast and can be reasonably accurate for certain types of molecules. CNDO/S and ZINDO/S were especially designed for p–p transitions in aromatic compounds, n–s, and other transitions may be obtained with a larger error. Finally, the TDDFT, for example, its linear response formulation,61 is based on the time-dependent KS equations 0 HKS ðtÞjCðtÞi ¼ i_ qqt jCðtÞi (cf. equation 4). The presence of a time-dependent term Vext(t) in the Hamiltonian causes a temporary density variation. Typically, TDDFT equations are formulated in a matrix form,62 such as 

A A

B B



o

C 0

0 C



X Q ¼ ; R Y

ð22Þ

    where Aia;jb ¼ dRij dab ðej  eb Þ=ðnb  nj Þ  Kia;jb, Bia,jb ¼  Kia,bj Cia;jb ¼ dij dab 1=ðnb  nj Þ , ei are orbital energies, the ‘coupling matrix’ Kia;jb ¼ dr dr 0  ji ðrÞja ðrÞfXC ðr; r 0 ; oÞjj ðr 0 Þjb ðr 0 Þ, fXC(r,r0,o) is the exchange-correlation kernel, ni are orbital occupations, and Q and R describe the perturbation. The corresponding TDHF method (for H0 KS ¼ H0 HF) is also called RPA (see Section 8.27.2.2). We can see that the order of matrices K, A, and B is a square of the product of the number of occupied and virtual orbitals, (occ  virt)2. For large molecules, the matrix elements must be computed on the fly. Equation 22 is solved iteratively. TDDFT is currently the most frequent approach for the computation of excited states, as it is reasonably accurate and fast. The TDDFT equations lead to diagonalization of large matrices, which is carried out using so-called Davidson iterative schemes. The excited states are found one after the other, starting from the lowest energy one. TDDFT ECD computation is frequently implemented in distributable software (e.g., Dalton, Gaussian, ADF, Turbomole). Hybrid functionals with HF exchange usually perform better than the GGA functionals. For example, B3LYP, PBE0, or BH and HLYP (with 20, 25, and 50% of HF exchange) may be recommended as the first choice. Also, recent functionals adapted for asymptotic behavior (CAM-B3LYP, wB97XD), double-hybrid functionals (B2PLYP), or functionals based on a statistical average of orbital potentials (SAOP) are expected to perform well for ECD. Similar to the ORD calculations, the environment should be involved, at least implicitly, using for example, PCM. The vibrational structure of absorption and ECD bands can be calculated most easily at the harmonic approximation.63 The vibrations are observed mainly on bands of aromatic compounds.

8.27.4.4

Origin Dependence of Theoretical ECD

Ab initio ECD spectra, if calculated in a limited AO basis set, can be dependent on the origin. Although a simple technique can be used to avoid it, the origin dependence is sometimes ignored, which can cause virtually unlimited computational errors. Therefore, it is important to at least schematically explain the problem. ^ ¼ ðe=2mÞr  p ^ ¼  ði_ e=2mÞr  r (or a We need to realize that the magnetic dipole operator in equation 17 is defined as m ^ 0 ¼ m ^ ði_ e=2mÞT  r. similar sum for more particles). Upon an origin shift by a translation vector T, the operator changes to m ^eT. For orthogonal states g and k, we obtain hgj^ ^ ¼  er changes to l ^0 ¼ l l0 jki ¼ hgj^ ljki and Similarly, the electric dipole operator l ^ 0 jgi ¼ hkjm ^ jgiði_ e=2mÞT  hkjrjgi, so that an extra term appears in the expression (17) for the rotational strength, hkjm R0 ¼ R þ ð_ e2 =2mÞhgjrjki  T  hkjrjgi.

530

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Obviously, the second term must be zero for exact wavefunctions, as the rotational strength is an observable quantity, independent of the arbitrary choice of origin in the calculations. The dipole-velocity transformation becomes important here: hgjrjki ¼

_2 hgjrjki; mðek  eg Þ

ð23Þ

which is valid for exact wavefunctions and energies ei. However, we can always calculate the dipole integral as a gradient (‘velocity’) using equation 23. Then, because of a vector identity A  B  A ¼ 0, the problematic term is zero by definition, also for approximate wavefunctions,

_ e2 _ 3 e2 hgjrjki  T  hkjrjgi ¼ hgjrjki  T  hkjrjgi ¼ 0: 2m 2ðek  eg Þm2 The velocity formalism should always be used when calculating ECD with field-independent AOs (e.g., in the Gaussian program). Another method of eliminating the origin dependence is TDDFT formulation using GIAOs, which is used, for example, in the Dalton software.

8.27.4.5

Applications

10

10 Ac-(Ala)4-NHMe

Ac-(Ala)4-NHMe Δ

Δ

5 0

−5 Ac-(Ala)6-NHMe

0

Ac-HAAAH-NH2 Ac-HAAAHAA-NH2

2

−10

0 Ac-(Ala)6-NHMe

−10 200

210

220  (nm)

230

240

4

200

210

220  (nm)

230

[Θ] (103 deg cm2 dmol−1)

A typical application of ECD consists of monitoring of conformational changes. But this technique found broader use in structural studies (determination of the absolute configuration, protein, and nucleic acid secondary structure), kinetic and thermodynamic studies of protein and DNA folding, aggregation studies (enantioselective interaction, protein–ligand, or DNA–ligand interactions, determination of the binding constants), etc. Despite the limits in accuracy, TDDFT simulations made the process significantly more reliable for a broad scale of organic compounds.24,64 The power of a combined experimental–theoretical analysis can be found in a BINOL (1,10 -binaphtyl-2,20 -diol) solution study.65 Several spectroscopic methods (ECD, VCD, OR, and IR) were used, and spectra were modeled on MD and DFT structures. For ECD, the B3LYP and CAM-B3LYP functionals were used; the performance of the latter was slightly superior. The absolute configuration of all three stereogenic axes and the molecular structure of 1,10 -binaphtyl-based phosphoramidites were determined. So far, ab initio and DFT methods have been too demanding to be used routinely for biopolymers. For proteins, experimental ECD spectra can be decomposed into curves of standard secondary structures.66 Many decomposition schemes were proposed, including factor analysis, singular value decomposition, and neural networks.67 The decomposition method is satisfactory for determining the a-helical content. Other forms (b-sheet, b-turn, 310-helix, random-coil structure) provide less distinct spectra.68 In our recent work, we explored several theoretical approaches to better understand various aspects of CD spectra of a-helical peptides,69 including the chain-length. The combined MD and quantum-chemical (TDDFT) techniques were able to capture the most important effects, such as the dynamics, hydration, and peptide length. However, a quantitative modeling remains a challenging task. TDDFT yielded similar dependence of the spectra on the peptide length as the simplified TDC model (see Figure 3). The modeling was consistent with experimental spectra of short peptides containing two histidine residues, where the a-helix was stabilized by the cobalt ion. The influence of the side chains and other protein parts on peptide ECD was predicted to be limited.

240 190 200 210 200 230 240 250  (nm)

Figure 3 Comparison of shorter and longer oligopeptide ECD spectra obtained by the TDDFT (left) and TDC (middle) methods with the experiment (right).69 Reproduced with permission from ACS.

In another DFT study, we modeled porphyrin ECD induced by a chiral environment.70 Several mechanisms by which ECD can be induced in a nonchiral chromophore (electrostatic perturbational, dipolar, and direct covalent) were discussed. TDC between

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

531

chirally stacked porphyrins was determined to be the most important contribution (cf. Figure 4). This is in agreement with previous experimental observations, where chiral matrices often induce the stacking and large ECD signals. Approximately 10 times smaller signal could be induced by a chiral orientation of the phenyl or similar residues covalently attached to the porphyrin core. Also, this prediction is in agreement with known experiments. Perturbation models realized by a chirally arranged porphyrin and a point charge or by a porphyrin and the methane molecule provided the smallest signals.

Δ ( l mol−1cm−1)

800

 d

60

600 45 67.5

400 200

22.5

90

0 −200 −400 300

250

350

400

450

 (nm)

d = 4.5 Å

1000

0 −1000 1000

d

d=6Å 0

Δ ( l mol−1 cm−1)

Δ ( l mol−1cm−1)

1000

d = 4.5 Å

0 −1000 1000

d=6Å

0 −1000

−1000 320

360 400  (nm)

440

320

360 400  (nm)

440

Figure 4 Orientation and distance ECD dependence for two porphyrin rings, as modeled by TDDFT.70 Reproduced with permission from ACS.

Another work deduces interaction types (ion pairing, hydrogen bonding or p–p stacking) of two dissimilar chiral molecules from a comparison of experimental and calculated ECD spectra.71 The authors investigated a tert-butylcarbamoylquinine selector and (S)-3,5-dinitrobenzoyl alanine selectand, using DFT methods for ECD of free selector and selectand as well for their complex. The complexation-induced changes could be modeled with a reasonably good agreement with experiment.

8.27.5

Magnetic Circular Dichroism as a Chiral Method for Achiral Molecules

MCD is a chiroptical method that uses circularly polarized light, but it does not need chiral molecules. Instead, a chiral environment is created by a static magnetic field and the direction of the light propagation. Most MCD applications involve highly symmetric systems (e.g., porphyrins), often with degenerate orbitals. The effect is measured as a difference in absorption of the left- and right-circularly polarized light exhibited by the sample in the presence of a static magnetic field oriented along the direction of light propagation. Although the effect occurs for electromagnetic radiation ranging from IR to X-ray, MCD most often refers to the UV–vis region. A vibrational variant using IR light is known as the magnetic vibrational circular dichroism (MVCD).72 X-ray MCD is predominantly used for solid-state studies.73 MCD spectroscopy can be conveniently carried out using an ECD instrument equipped with a permanent magnet, an electromagnet (both with fields up to 1.5 T), or a superconducting magnet (up to B8 T). Thus, many suppliers of ECD spectrometers also sell the magnets. Note that both MCD and ECD are measured as the absorption differences and both involve electronic transitions. Usually, the experimental conditions are also similar.

532

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

8.27.5.1

Theory of MCD

For a detailed derivation of the expressions and their discussion, we refer the reader to Reference 74. The CD absorption difference associated with a transition from a ground state g to an excited state j is De ¼ l

2  2 Ng  Nj  g9^ mþ 9j    g9^ m 9j  rðEÞ; N

ð24Þ

where Ng and Nj are numbers of molecules in the ground and excitedpstate, respectively, and N is the total number of molecules in ffiffiffi ^ is the electric dipole moment operator. The light is ^7 ¼ ð^ unit volume, l is a constant independent of energy, m mx 7^ my Þ= 2 and l propagating along the z-axis. Equation 24 can be obtained from the Fermi golden rule in perturbation calculus assuming the dipole approximation. Apart from the electromagnetic perturbation, we need to introduce the perturbation by the static magnetic field B. According to the time-independent perturbation theory (cf. equation 9), we express states g and j with the aid of field-free wavefunctions g(0) and k(0), for example,      ð0Þ E X kð0Þ mB ^ gð0Þ  ð0Þ E  k  ; ð25Þ jg i ¼ g ð0Þ ð0Þ Eg  Ek ka g ^ is the magnetic dipole moment operator and E(0) where m are the unperturbed energies. i We assume that the band shape provided by the density of states r(E) does not change in the presence of the magnetic field, but only shifts in energy, which is referred to as the rigid shift approximation. This shift is a function of the magnetic field, so that the density of states in the presence of a magnetic field is ^ jgiB þ hjjm ^ jjiB : ð26Þ rðEÞ ¼ r Eð0Þ  hgjm Finally, we need to take into account possible splitting of degenerate ground-state energy levels and their Boltzmann probabilities pg. For a weak magnetic field, ^ jgiB=ðkTÞ: pg B exp½Eg =ðkTÞE1 þ hgjm

ð27Þ

Combining this knowledge, we can write the MCD absorption difference as a sum of three parameters, so-called Faraday A, B, and C terms.42,74 For an isotropic sample (solution), we obtain



qr 1 r; ð28Þ þ BþC DeB A  kT qE where A¼

eabg X D ð0Þ   ð0Þ E D ð0Þ   ð0Þ E D ð0Þ   ð0Þ ED ð0Þ   ð0Þ E ^a j ^a g ^b j ^g g  g m Im g m j m ; j m 2n g "    ð0Þ  D X kð0Þ m   ð0Þ ED ð0Þ   ð0Þ E ^ a g eabg X ^b j ^g k B¼ j m Im gð0Þ m Ek  Eg n g  ð0Þ ka g ð0Þ  D X j 9m   ð0Þ ED ð0Þ   ð0Þ E ^ a 9k ^b j ^g g k m ; gð0Þ m þ Ek  Ej ka j eabg XD ð0Þ   ð0Þ E D ð0Þ   ð0Þ ED ð0Þ   ð0Þ E ^a g ^b j ^g g C¼ Im g m j m ; g m 2n g

and n is the degeneracy of the ground state g. We can see that the A term is associated with the degeneracy of the ground and excited states (otherwise, they would not have a permanent magnetic moment), that is, it is zero when neither a ground nor an excited state is degenerate. The C term also needs a degenerate ground state, but originates in different populations of the split ground-state electronic levels. The B term originates from the mixing of the states induced by the magnetic perturbation. Quite often, the C term is neglected, and MCD is considered as a sign of the degenerate states entering term A, which leads to a conservative split signal (proportional to q r/q E, see equation 28) or a transition affected by mixing with other states, which leads to one-sign signal B (proportional to r(E)). It is the B term that enables us to study virtually any molecule by MCD because the effect of the mixing of states by a magnetic field is present to some extent in all matter. However, one has to realize that this classification of MCD terms is somewhat arbitrary and cannot be applied for nearly degenerate cases. The theory must be modified, for example, in case of breaking the degeneracy by spin-orbit coupling.75 Apart from the quantum chemical calculations, predominantly using TDDFT, MCD spectra are sometimes analyzed using semiempirical approaches. For example, a simplified MO method was proposed for interpreting MCD spectra of porphyrin-based

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

533

molecules.76 It was found useful for understanding the structure–spectrum relationships; obviously, its applicability outside of the field of the cyclic aromatic compounds is limited.

8.27.5.2

Practical Calculations

Because of the complexity of the MCD phenomenon, it is not yet routinely calculated, and the choice of suitable software is limited. Currently, it is available in the Dalton and ADF programs. Dalton enables to calculate either a pure B term, using the keyword ‘MCDBTERM,’ or a combined contribution of the A and B terms, utilizing the polarization propagator approach.77 The latter method includes a damping parameter needed in the input. This enables a direct simulation of the band shape. Both methods are implemented at the TDDFT level. The ADF program supports calculations of separate A, B, and C terms, also using TDDFT. Note that currently, ADF cannot use hybrid functionals for MCD. In Figure 5 we document how a porphyrin protonation can be monitored by MCD, and reproduced by the TDDFT computation implemented in ADF (unpublished results). The higher wavelength (B650 nm) part is calculated with a larger error than the lower one, which is typical for such computations. Otherwise, the main changes in the intensity and wavelength observed experimentally with lowering of pH are well reproduced by the theory. pH = 2.5 (H2TPPS2−)

Exp.

ΔA-MCD

Calc.

Calc.

pH = 6 (TPPS4−) Exp. 300

400

500

600  (nm)

700

800

TPPS4−

H2TPPS2−

Figure 5 Experimental MCD spectra of tetraphenylporphyrinsulfate (TPPS4, and its protonated form H2TPPS2), and SAOP/TZP/COSMO calculations for TPP.

8.27.6 8.27.6.1

Vibrational Circular Dichroism Experimental Aspects

Optical activity caused by vibrational motions was applied for structural studies later than the chirality caused by electronic transitions (ORD, MCD, and ECD). One of the main reasons was the weakness of the VOA phenomena, given by the relative importance of the magnetic and electric dipolar terms in the expansion of the interaction of molecules with the electromagnetic field. Although ORD in infrared had been known for a long time,78 VCD was measured much later, in 1973, on inorganic crystals,79 and in 1974 on 2,2,2-trifluoro-1-phenylethanol.80 Typically, the ratio of VCD to absorption (equation 20) is approximately 10 times weaker than for ECD. However, thanks to advances in the instrumentation (lock-in amplifier, dual polarization modulation, careful calibration)81 the technique could soon be applied to a large number of chiral molecules in the liquid state, but also solid suspensions and gases (see Chapter 8.25). The driving force of the VCD development was the richer information of the IR spectra than for UV–vis. Although dissolved molecules have often only one or a few broad absorption bands in the UV–vis region, the vibrational transitions are more numerous (3N–6 fundamental modes for a molecule with N atoms) and better resolved. The longer wavelength in VCD has some disadvantages (lower dissymmetry g ratio than for ECD), but also enables measurements on macromolecular systems that are not well suited for ECD because of the light scattering.82 Finally, the computation of vibrational properties is easier than that for the electronic ones in that it does not require electronic excited states.

8.27.6.2

VCD Theory

The experimental results required a rigorous interpretation. Originally, simplified approaches were developed, such as the fixed partial charge (FPC) model of Schellman.83 To calculate IR absorption and the VCD intensity of a vibrational transition90S-9iS,

534

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

we need the dipole and rotational strengths, the same ones as in ECD,  Di ¼ li  li ¼ h0j^ ljiij2 ; ^ j0i; Ri ¼ Imli  mi ¼ Imh0j^ ljiihijm

ð29Þ ð30Þ

P ^ are the respective electric and magnetic dipole operators. Within FPC, we write the moments as l ^ and m ^ ¼ l0 þ A qA DRA where l P ^ ¼ m0 þ 1=2 A qA RA  PA =mA , where qA are the atomic partial charges, mA, RA, and PA are nuclear masses, positions, and and m momenta, and DRA is the vibrational deviation of atom A from the equilibrium position. Without describing the FPC model further, we can say that it is no longer used. However, it was useful to gain an idea of the size and origin of molecular VCD. It provided the correct magnitude of the phenomenon,84 but failed in many cases to reproduce the sign of VCD bands, which is the most important information provided by the spectra. The failure of the model can partially be attributed to the lack of a method to generate accurate vibrational frequencies in the origins of VCD (B1970). More importantly, however, the concept of the fixed charges proved to be physically incorrect. Electronic currents strengthening or weakening the movement of the nuclear charge accompany molecular vibrations, which cannot be predicted without proper quantum chemical computations. The more general coupled oscillator (CO, see also Section 8.27.4.2) model85 was more successful, and is still being used in approximate computations. The model is simple, exact for separated chromophores,86 and amenable to many systematic improvements.49 Let us recall that as input parameters we need transition electric dipole moments li located on chromophores at positions ri and transition frequencies oi. Then we set out the Hamiltonian interaction matrix H, where the diagonal elements are Hii ¼ _ oi and the interaction is approximated by the dipole–dipole potential, Hij; ja i ¼ ð4pe0 Þ1 rij2 li  lj  3rij  li rij  lj =rij5 . After diagonalization of H, we obtain the transition energies eJ and eigenvectors cJ, so that the resultant dipole and rotational strengths will be (in atomic units) 2 X   ð31Þ DJ ¼  i cJi li  ; and RJ ¼ 12

X X cJi li  ej cJj rj  lj : i

ð32Þ

j

P The model predicts conservative spectra, J RJ E0, such as the couplets for the C ¼ O stretching modes observed in a-helical proteins or nucleic acids. For more rigorous simulations, we need quantum mechanics. Within the harmonic approximation a VCD computation can be thought of as having two major parts: computation of the force field (FF, from which we obtain the transition frequencies) and spectral intensities. Therefore, we should briefly review the properties of the harmonic vibrational Hamiltonian, ! 3N ^ 2 3N X 3N X Pi 1 X 1 t ^ ^ ^ H¼ þ Fij DRi DRj ¼ ðp p þ qt  f  qÞ; ð33Þ 2 i ¼ 1 mi i ¼ 1 j ¼ 1 2 ^ i , masses mi, the force field Fij, also referred to as the Hessian, and deviations from equilibrium which consists of nuclear momenta P pffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi nuclear positions DRi. We wrote it also using mass-weighted coordinates qi ¼ mi DRi and momenta ^pi and FF fij ¼ Fij = mi mj. In the early days (mid twentieth century), empirical force fields were proposed to explain the spectra, often refined according to the experiment. Today, most quantum chemical software packages provide the FF computed as second energy derivatives (Section 8.27.2, Table 1), either numerically from gradients or energies or more accurately and faster using the CP techniques.87 ^ k ), To simplify the Hamiltonian, we introduce so-called vibrational normal t mode coordinates (Qk) and associated momenta (P P P 1 ^ k ¼ 3N s1 ^p , by which we obtain H ^  st  s  P ^ þ Qt  st  f  s  Q . When the transformation is ^ s q and P ¼ 1=2 P Qk ¼ 3N j j¼1 kj j¼1 kj j unitary (st  s ¼ 1, where 1 is the identity matrix) and the transformed force field, st  f  s ¼ K, is diagonal (i.e., Lij ¼ 0 for iaj and potential the Hamiltonian becomes a sum of one-dimensional Lii ¼ o2i ), which can be always done for a quadratic form, P ^ P 2 ^ þ o2 Q2 ¼ 3N h ^ ¼ 1=2 3N P ¨harmonic oscillator Hamiltonians H i i i i¼1 i¼1 i ðQi Þ. The one-dimensional harmonic oscillator Schro dinger equations

1 q2 2 2 _ 2 þ o Q ð34Þ i i wi ðQi Þ ¼ Ei wi ðQi Þ 2 q Q2i have analytical solutions with energies Ei,vi ¼ _ oi(1/2 þ vi), quantum numbers vi ¼ 0,1,2,y, and wavefunctions wi. Most often, we are only interested in so-called fundamental transitions, from a vibrational ground state (all vi are zero) to a first excited state (for some j vj ¼ 1). Because of the translations and rotations, there are in general 3N–5 (for linear) and 3N–6 (for other molecules) vibrational modes with nonzero energies. Note that these were actually obtained without solving the differential Schro¨dinger equation, only by coordinate transformation and diagonalization. Instead of s, a direct Cartesian–normal mode transformation matrix S can be defined as Skj ¼

q DRj pffiffiffiffiffi ¼ skj = mj : q Qk

ð35Þ

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

535

To obtain the IR absorption intensities, it is convenient to expand the electric molecular dipole in the normal mode dis P placements ma ¼ ma ð0Þ þ j q ma =q Qj 0 Qj þ y, so that the transition dipole element for a fundamental transition j is sffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffi

_ q ma _ X A ¼ P SAb;j : ð36Þ hw1 jma jw0 i ¼ 2oj q Qj 0 2oj A;b ba   A ¼ q ma =q RAb 0 ¼ EAba þ ZA edba are also called the ‘atomic polar tensors‘ (APT, according to the The ‘dipole derivatives’ Pba transformation properties under coordinate inversion) and have electronic and trivial nuclear parts; ZA is the charge of nucleus A. Similar to the FF, APT is available in quantum chemical programs, most conveniently calculated using the CP approaches, for example, as *   +   qj

  g A ^  Eba ¼ 2 jg m ; ð37Þ  a  q RAb 0 ^ a. that is, from a ground-state wavefunction jg, its derivative, and the electronic dipole moment operator m For the rotational strength needed for the VCD intensities, we cannot expand the magnetic moment in the coordinates, as q ma/ q Qj would be zero in the BO approximation. However, we can consider expansion  in nuclear momenta, which are formally P equivalent to coordinates in the Hamilton formalism, ma ¼ ma ð0Þ þ j q ma =q Pj 0 Pj þ ?. Indeed, this does lead to the right contribution; within the magnetic field perturbation (MFP) theory,6 the transition magnetic dipole moment can even be expressed using only the electronic ground state. For a fundamental transition of mode j,  1=2 X A Mba SAb;j ; ð38Þ hw1 jma jw0 i ¼  2_ 3 oj A;b

where MAba ¼ (iMA/2_)(q ma/q PAb)0 ¼ IAab þ (ieZA/4_)eagbR0Ag is the atomic axial tensor (AAT), R0A is the equilibrium position. Note that it can be considered as a magnetic dipole derivative, but with respect to the momentum, not a coordinate. The electronic part can be obtained within the CP-MFP formalism as 



 q jg  q jg A  : ð39Þ Iba ¼  q RAb 0 q Ba 0 A response–function formalism has emerged recently, suitable for more efficient computer implementations.88 Alternative approaches, such as the vibrational coupling theory (VCT)89 or an excitation (SOS) model,90 are less accurate and are not used nowadays. The original MFP formulation had the limitation of origin dependence of the resultant rotational strength if calculated using approximate wavefunctions (limited basis set of atomic orbitals). This could be circumvented by the distributed origin gauge (DOG) transformation;91 however, computation with GIAOs is more universal and more frequent (cf. Section 8.27.2.2).92

8.27.6.3

Spectral Simulations

The simulations of VCD using the harmonic approximation and the CP-MFP-GIAO method have become standard within DFT, and provide satisfactory results for medium-sized molecules. The expansion of the methodology beyond the harmonic approximation is possible;93 however, the anharmonic contribution within the most experimentally accessed range (B1100–1800 cm1) is usually small compared with other experimental and theoretical uncertainties. Modern functionals provide reasonably accurate frequencies and spectral intensities; a typical example is shown in Figure 6, where the prevalent conformer of a relatively complicated tetramer of S-2,20 -dimethyl-biphenyl-6,60 -dicarboxylic acid was determined using VCD spectroscopy.94 Hybrid functionals (B3LYP, B3PW91, etc.) seem to be slightly more reliable than GGA ones (BP86, BPW91).95 Accurate quantum-mechanical VCD simulations of larger molecules are not possible. However, in this case, locality of the vibrational phenomena can be explored. For example, on a simplified peptide dimer, we could reproduce the main VCD spectral features of standard peptide secondary structures.96 The coupling between individual chromophores important for VCD is primarily mediated by the spring-like short covalent bond interactions. Recently, the idea was generalized in the Cartesian coordinate transfer (CCT)97 of molecular property tensors, such as FF, APT, and AAT. The tensor elements, for example, FF dependent only on two atoms, are calculated for a smaller fragment and then transferred onto a target molecule. The concept of force field transferability has been known for a long time. However, the CCT implementation proved to be more convenient for automatic computer manipulations than the classical approaches using manually defined internal coordinates. Also, within CCT, the properties of the large molecule can be iteratively improved up to the exact values. A very useful enhancement of the method was the partial normal mode optimization,98 which allows one to relax spectroscopically important coordinates, but does not introduce unnecessary geometry deformation. In this way, spectra of large proteins and nucleic acids were calculated, as exemplified in Figure 7 for a model (dG-dC)4 octamer.

536

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Δ a

0.05

−0.05 a Expt Calc

b



500 a

1700

1600

1500

1400

1300

1200

1100

Figure 6 Example of an agreement between a CP-MFP-GIAO (B3LYP/6-31G) calculation and experiment, for VCD (De) and absorption (e) spectra of a tetramer complex.94 Reproduced with permission from ACS.

Current challenges of the computational VCD methodology comprise the proper involvement of the solvent and reasonable averaging of many conformational minima in flexible molecules. For strongly polar solvents generating hydrogen bonds, the usual PCM dielectric continuum approximation is not adequate.100 The flexibility can be taken into account, for example, by a combination of classical and quantum molecular dynamics.101 The molecular motion may be observable by VCD in the future more directly by time-dependent experiments; new cross-polarization detection schemes enable the simultaneous measurement of VCD and ORD in the vibrational region.7 The magnetic alternative of VCD102 is not in use much as yet, neither is the rotationally resolved magnetic VCD for gases.103

8.27.7

Raman Optical Activity

8.27.7.1

History and Experiment

The ROA is a difference in scattering of the right- and left-circularly polarized light, IROA ¼ IR  IL. We note that for historical reasons, the ROA intensity is defined in the opposite sense (‘right-’ minus ‘left-’) than for the dichroic (CD, VCD) differential methods. The first molecular ROA was measured in 1973,104 only slightly earlier than VCD. But these two forms of VOA exhibit significant differences. Experimentally, ROA is easier to measure for aqueous solutions and in a broader range of wavenumbers than VCD. This makes it a convenient tool for many biologically relevant systems. The laser light (typically at 532 nm) used for the Raman excitation is suitable for conventional optical elements. However, ROA samples need to be free of impurities that cause fluorescence and inhibit the measurement, but which would not interfere with the IR (VCD) experiments. Various experimental ROA setups are possible, with forward, right angle, and backward scattering, with a circular polarization of the incoming light, a dual polarization, etc. (see Chapter 8.26).105

8.27.7.2

ROA Theory

Theoretically, ROA is a two-photon scattering process, unlike the one photon VCD. Higher order frequency-dependent molecular polarizabilities are responsible for ROA of an isotropic sample. A dipole is induced in a molecule in the electromagnetic field, 1 1 q Bi;b þ ?; mi;d ¼ ai;da Ei;a þ Aid;bg rEi;bg þ G0idb 3 o qt

ð40Þ

537

1690

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

1710 1678

Calculation

B-form

Experiment

1691

(b)

1707

Δ

1693

Calculation

1671

Z-form

1656

Experiment

1600 (a)

1400

1200

Wavenumber

1000

800

(cm−1)

(c)

Figure 7 Comparison of the calculated and experimental absorption and VCD spectra (a) for the B- (b) and Z-(c) forms of DNA, for a model (dG-dC)4 octamer. The CCT calculation99 reproduced the experimentally observed sign flip of the couplet at B1700 cm1.

where E and B are the electric and magnetic fields, and a, A, and G0 are the respective electric dipole–electric dipole, electric dipole–electric quadrupole, and electric dipole–magnetic dipole polarizabilities. Note that from the G0 tensor, we could obtain the OR (Section 8.27.3). G0 contains the product of the magnetic and electric dipole, as in the rotational strength; thus, the CD theory can also be reformulated using this tensor. For nonresonance ROA, we can use the Born–Oppenheimer approximation (also referred to as Plazcek’s approximation in this context). The tensors for a molecular electronic ground state 9gS are expressed as sums over electronic excited states 9kS, X hgj^ ma jkihkj^ mb jgi aab ¼ 2 okg Re ; 2  o2 o kg ka g G0ab ¼ 2o

X hgj^ ^ b jgi ma jkihkjm Im ; 2  o2 o kg ka g

ð41Þ

ð42Þ

and _

X ma jkihkjY bg jgi hgj^ Aa;bg ¼ 2 okg Re : o2kg  o2 ka g

ð43Þ

_ P The electric quadrupole can be defined (the definition varies slightly in different sources) as Y ¼ 1=2 i qi ðri ri  1ri2 =3Þ. Furthermore, we need to introduce derivatives of these tensors with respect to the normal mode coordinates, for example, 0i =9, ai ¼ qa/qQi, and so-called isotropic tensors invariants, a2i ¼ aiaa aibb =9, b2i ¼ ½3aiab aiab  aiaa aibb =2, aG0i ¼ aiaa Gbb 0 i 0i i 0i i i bGi ¼ ½3aab Gab  aaa Gbb =2, and bAi ¼ oaab eagd Ag;db =2. We use the Einstein summation convention, that is, indices occurring twice in a product are summed over. For isotropic samples (solutions) within the harmonic approximation, we can then express

538

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

the Raman (IRam ¼ IR þ IL) and ROA (IROA ¼ IR  IL) intensities for different experimental setups.42,105,106 As an example, for the backscattering (1801) experiment, we obtain IR þ IL B 90a2i þ 14b2i ;

ð44Þ

IR  IL B ð2=cÞ½24bG0i þ 8bAi :

ð45Þ

Note that the absolute intensity is not measured for scattering, unlike for VCD/IR. However, the ROA/Raman circular intensity difference ratio (CID ¼ IROA/IRam) is an important experimental parameter that can also be compared with the calculation. The vibrational frequencies are obtained in the same way as for VCD; the normal mode tensor derivatives are usually not calculated directly, but obtained from the Cartesian derivatives and the transformation matrix (equation 35). For simulation of the theoretical spectra, it is important to consider the Boltzmann populations of vibrational states. The Raman scattering techniques allow monitoring of vibrations with very low frequencies (comparable to the temperature quantum, B200 cm1) involving transitions from excited states that are more intense than those from the ground state. We can also use the Lorentz shape (with a width D) to simulate inhomogeneous band broadening; thus, the calculated Stokes spectrum to be compared with experiment is 1 h o i1 1  o  o 2 i i 4 þ1 ; SðoÞ ¼ IRam=ROA 1  exp  oi kT D

ð46Þ

where k is the gas constant and T is the temperature. Note that neither the Boltzmann factor nor the division by the vibrational frequency (1/oi, coming from the normal mode transition integral) is included in the intensities provided by common quantum chemical programs, such as Dalton,107 Gaussian,108 Turbomole,109 and CADPAC.21 Empirical schemes for computation of the ROA intensities are relatively rare, such as Barron’s two-group model42 and its extension in the polarization model,110 which uses the dependence of the G0 and A tensor derivatives on the origin:

q Gab q Gab io q aad ebgd Rlg ð0Þ ¼ ðlÞ  ; ð47Þ 2 q Rle q Rle q Rle   q Aa;bg q Aa;bg 3 l q aag q aad l q aab R ð0Þ ¼ ðlÞ þ þ R :  dbg Rld b g 2 q Rle q Rle q Rle q Rle q Rle

ð48Þ

The irreducible contributions expressed with the origin at the differentiated atom l (q Gab/q Rle (l) and q Aa,bg/q Rle (l)) are neglected, and ROA is thus obtained from the common origin (0) representations using only the a-tensor. The contribution of the polarizability a to ROA is usually dominant, and the model allowed for a fast estimation of the spectra, for example, for a very large set of conformers.111 After G0 and A became available with the same computational effort as a, the polarization model became obsolete. The first full quantum chemical ROA computations implemented in the CADPAC program at the CPHF level involved a P relatively minor simplification, the static limit of the G0 -tensor. For o-0, we obtain G0ab =o ¼ 2 Ea G o2 EG  E B ^ b jGiB Ca Cb S. This means that the tensor could be calculated as an overlap integral of electric (E) and magnetic ImhGj^ ma jEihEjm (B) field derivative of the ground-state electronic wavefunction. The origin dependence of the results can be resolved by a dipole-velocity transformation (equation 23) for the a-tensor, similar to the DOG gauge in VCD. Soon, however, the origin-independent ROA formulation in GIAOs emerged.112 Today, frequencydependent analytical ROA tensor derivatives can be obtained, for example, within the response theory.106 Note that the variables needed to calculate VOA (ROA and VCD) are also summarized in Table 1. For larger molecules, the CCT scheme97 can be used; a special fragment approach was proposed for peptides.113 A normal mode tracking technique and selection of intensity carrying modes114 can also save computational time if only a limited part of the spectrum is required.

8.27.7.3

Simulation Examples

ROA experiments and computations, similar to VCD, are useful tools for the determination of absolute configurations for small and rigid molecules. In this case, a very good agreement can be achieved between the simulated band intensities and frequencies. Thus, for example, the absolute configuration of bromochlorofluoromethane115 could be determined by ROA, probably as the only possible technique. The harmonic approximation is usually sufficient to provide unambiguous band assignment in the usual interval of wavelengths (B200–2000 cm1), although anharmonic corrections may be useful.95 This is documented in Figure 8 for S-methyloxirane; note that the vibrational configuration interaction (VCI) anharmonic correction clearly improved both simulated frequencies and intensities for wavenumbers 4B1200 cm1.

IR − IL

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

4 2 0 −2 −4 4 2 0 −2 −4

539

Harmonic

Anharmonic

Experiment

400

600

800

1000

1200

1400

1600

Wavenumber (cm−1) Figure 8 Example of CP-DFT (B3LYP/aug-cc-pVTZ) ROA computation for S-methyloxirane, in the harmonic and anharmonic approximation, and experimental spectrum of a neat liquid.

For flexible molecules, a conformer averaging is needed, both for well-defined conformers116 and for flexible molecules oscillating around an equilibrium.117 Molecular flexibility in general makes the ROA signal smaller and causes a large inhomogeneous broadening of Raman and ROA spectral bands. As an advanced example, we can see the agreement between the calculated and the experimental spectrum of the b domain of rat metallothionein in Figure 9.118 The computation well captures the most intense spectral features, such as the couplet at B1300 cm1. Above 1400 cm1, however, the precision of the computation is limited, most probably because of the anharmonic or solvent effects, or fine structural inaccuracies, and does not agree with the experiment.

(I R − I L) (arb. units)

(I R − I L) (A4 a.m.u.−1)

0.004

800

1000 1200 1400 1600 Wavenumber (cm−1)

0.003 0.002 0.001 0 −0.001 −0.002 −0.003

800 1000 1200 1400 1600 Wavenumber (cm−1)

Figure 9 Experimental (left) and calculated (BP86/RI/TZVP) ROA spectrum of the b domain of rat metallothionein.118 Reproduced with permission from ACS.

For the potassium–valinomycin complex (Figure 10) computed spectra of more conformers could be compared with the experiment. In this case, a combined direct DFT-CCT strategy was chosen for the calculations, since limited accuracy was found for transferred ROA tensors. The conformer ratios, similar to model dipeptides,120 agreed with those obtained by NMR spectroscopy. Some conformers cannot be resolved by NMR due to a fast exchange; in such cases, the ROA methodology is a unique tool for structural studies in solutions.

8.27.8

Common Aspects of the Computational Approaches

We saw that the chiroptical spectra could be analyzed efficiently with the use of the quantum chemical methods. The most common computational programs providing molecular optical activity are summarized in Table 4. In Table 5, other state-of-theart software projects are indicated, although we realize that the software list cannot be complete, at least due to the dynamic development in the field.

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

( l R + l L ) / 1010

( l R − l L ) / 106

540

5 0 −5

Exp. Calc.

2 1 0 1800

1400

1000

600

200

Wavenumber / cm−1 Figure 10 Comparison of experimental and calculated Raman and ROA spectra of a potassium valinomycin complex. For the computation, B3LYP/6-311 þ þ G/PCM FF was obtained from six overlapping fragments.119 Reproduced with permission from the PCCP Owner Societies.

Table 4

Selected computer programs providing chiroptical properties

Program

Present capabilities

121

ADF www.scm.com CADPAC21 www-theor.ch.cam.ac.uk/software/cadpac.html Dalton,107 dirac.chem.sdu.dk/daltonprogram.org Gaussian,108 www.gaussian.com/ Turbomole109 www.turbomole.com/

Table 5

MCD, VCD, CD, ORD, ROA oriented to GGA functionals, STO basis sets VCD, ROA ROA, CD, VCD, MCD, ORD large flexibility in computational parameters, difficult to learn ROA, VCD, CD, ORD ‘the standard’ in computational chemistry CD, local versions also ROA118 efficient code, limited number of molecular properties

Other common ab initio software

Program

Link

ACES/CFOUR

http://www.qtp.ufl.edu/Aces2/ http://www.cfour.de/ http://www.msg.chem.iastate.edu/gamess/ http://www.cfs.dl.ac.uk/gamess-uk/ http://www.schrodinger.com/ http://www.teokem.lu.se/molcas/ http://www.molpro.net/ http://www.nwchem-sw.org/ http://www.thch.uni-bonn.de/tc/orca/ http://www.q-chem.com/

GAMESS Jaguar MOLCAS MOLPRO NWCHEM ORCA Q-Chem

Once formulated as a problem within the quantum theory, the OA simulations follow the general status of the computational apparatus development. The choice of the basis set, for example, is a minor problem in terms of the formalism, but very important for practical computations. Accurate results can be obtained only in the limit of an infinite basis set. Except for the ADF program, which uses Slater-type atomic orbitals, most programs exploit Gaussian-type functions. Plane wave and combination approaches could not be covered in this review, but may be a promising alternative for condensed-phase simulations. The theory can be

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

541

adapted for heavy elements where relativistic corrections are important, using approaches based on the Dirac equation. For example, core pseudopotentials (ECP, ‘effective core potentials’) can be introduced and used within a nonrelativistic apparatus. The accuracy ladder is well defined for the wavefunction approaches; following approximately the sequence HF, CIS, MP2, MRSCF, CCSD, CCSD(T)yFCI, we can systematically improve the accuracy. In particular, the CCSD(T) method is considered as the most versatile computational method providing experimental accuracy for most molecular properties. As pointed out above, the FCI results in an infinite basis are identical to the exact solution of the Schro¨dinger equation. Unfortunately, these methods scale sharply with the size of the system; formally, the time needed to calculate HF energies is proportional to the fourth power of the number of atoms (BN4), whereas the FCI even exhibits an exponential increase (BeN). Note that computations of properties are more demanding than those of energy only. For example, FF (second energy derivatives) require first derivatives of wavefunction, and scale like N5 for the HF method. The DFT methods are more convenient, typically with N3 (GGA) or N4 (hybrid GGA-HF functionals) scaling, but they are less reliable and less predictable than wavefunction computations. Therefore, DFT functionals need to be carefully tested before applying to a new problem. For example, functionals providing good vibrational frequencies (such as BPW91) may fail for ECD (where B3LYP is more reliable). A functional improvement for ECD (such as CAM-B3LYP) may deteriorate the performance for vibrations. The inclusion of the solvent in quantum chemical computations became possible only relatively recently, due to improvements in computational technology and computer speed. The solvent is an excellent example of a presumably minor aspect, neglected in the early days (1970–80s) of chemical computations, which became an annoying bottleneck as the accuracy standards increased. This is particularly true for chiroptical properties, sensitive also to the fine structure of hydration shells.

Initial geometry estimation X-ray, NMR structure, etc.

Energy calculation, best conformer selection Systematic conformer search, MD methods for larger molecules

Choice of the computational level for properties Important: for vibrations calculate force field at the same level as geometry

Tuning of computational parmeters

CCSD(T), augcc-pVTZ, HF, MP2, MP4, 6-311++G**, BPW91, B3LYP

Spectra simulation Including solvent and dynamic effects, anharmonicities...

Unsatisfactory agreement Comparison with experiment

Figure 11 ‘Daily life’ of a computational spectroscopist.

Good agreement

542

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

We can perhaps expect a similar need for proper accounting of molecular flexibility, dynamics, and anharmonic interactions. Fortunately, the solvent ‘ladder’ is relatively straightforward, gradually improving the ‘single-molecule-in-vacuum’ approximation by the point-charges model, the Onsager dipole model, COSMO and PCM dielectric continuum corrections, and, finally, clusters of explicit solvent molecules treated either as an MM part in the QM/MM approach or even as a part of the quantum system. Somewhat more problematic is the inclusion of the dispersion (van der Waals) forces, not very well described by HF or standard DFT, but extremely important for a correct description of many systems studied by OA (e.g., base pair stacking in DNA). Although advanced dispersion-corrected DFT schemes have became available in the past decade,122 we still find it difficult in practical computations to balance the dispersion with the solvent and dynamical aspects. Typical steps taken when simulating molecular chiroptical properties resemble those needed for other computations. Perhaps only the conformational aspect is more important, as the spectra are very sensitive to molecular shape. Yet, it may be useful to summarize the standard procedure (Figure 11), as a systematic and in the beginning a seemingly more laborious approach often pays off and yields correct results faster than ad hoc computations.

8.27.9

Conclusions

The interpretation of chiroptical spectra from first principles is possible due to the increased accuracy of the quantum chemical methods. Implementation of particular molecular properties in the quantum chemical software packages is driven by the advances in experimental instrumentation. As a feedback, reliable theory is encouraging the experimental methodology. Most common computations covered in this chapter use the density functional theory, which yields usable results for the interpretation of ORD, ECD, MCD, VCD, and ROA. The theory thus provides the link between the spectra and molecular geometry. The structure determination is complicated by the limited precision of theory, and approximate models of the solvent and dynamics. However, the complex theoretical descriptions required in spectral simulations provide precious information about the studied systems, not immediately apparent from experiment. The remarkable progress of computational chemistry and chiroptical spectroscopy in the last decades indicates that the methodology will remain a valuable tool for chemists and structural biologists interested in molecular structure, flexibility, and affinity to the environment.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33.

Harris, R. A.; Jameson, C. J. J. Chem. Phys. 2006, 124, 096101. Pasteur, L. The`ses de chimie et de physique; Bachelier: Paris, 1847. Craig, D. P.; Thirunamachandran, T. Molecular Quantum Electrodynamics; Dover: New York, 1998. Lee, T. D.; Yang, C. N. Phys. Rev. 1956, 104, 254–258. Grunenberg, J. Computational Spectroscopy; Wiley: Weinheim, 2010. Stephens, P. J. J. Phys. Chem. 1985, 89, 748–752. Rhee, H.; June, Y. G.; Lee, J. S.; et al. Nature 2009, 458, 310–313. Andrushchenko, V.; Bourˇ, P. J. Comput. Chem. 2008, 29, 2693–2703. Kohn, W.; Becke, A. D.; Parr, R. G. J. Phys. Chem. 1996, 100, 12974–12980. Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133–A1138. Stevens, R. M.; Lipscomb, W. N.; Pitzer, R. M. J. Chem. Phys. 1963, 38, 550. Olsen, J.; Jorgensen, P. J. Chem. Phys. 1985, 82, 3235–3264. McLachlan, A. D.; Ball, M. A. Rev. Mod. Phys. 1964, 36, 844–855. Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997–1000. Ditchfield, R. Mol. Phys. 1974, 27, 789–807; Wolinski, K.; Hinton, J. F.; Pulay, P. J. Am. Chem. Soc. 1990, 112, 8251–8260. Polavarapu, P. L.; Chakraborty, D. K.; Ruud, K. Chem. Phys. Lett. 2000, 319, 595–600; Kondru, R. K.; Wipf, P.; Beratan, D. N. J. Phys. Chem. A 1999, 103, 6603–6611. Pecul, M.; Ruud, K.; Rizzo, A.; Helgaker, T. J. Phys. Chem. A. 2004, 108, 4269–4276. Marchesan, D.; Coriani, S.; Forzato, C.; et al. J. Phys. Chem. A 2005, 109, 1449–1453. Rosenfeld, L. Z. Phys. 1928, 52, 161–174. Amos, R. D. Chem. Phys. Lett. 1982, 87, 23–26. Amos, R. D. CADPAC: The Cambridge Analytic Derivative Package, Issue 6.0; SERC Laboratory: Daresbury, UK, 1995. Costante, J.; Hecht, L.; Polavarapu, P. L.; Collet, A.; Barron, L. D. Angew. Chem. Int. Ed. 1997, 36, 885–887; Bak, K. L.; Jorgensen, P.; Helgaker, T.; Ruud, K. Faraday Discuss. 1994, 99, 121–129. Cheeseman, J. R.; Frisch, M. J.; Devlin, F. J.; Stephens, P. J. J. Phys. Chem. A 2000, 104, 1039–1046. Grimme, S.; Furche, F.; Ahlrichs, R. Chem. Phys. Lett. 2002, 361, 321–328. Stephens, P. J.; Devlin, F. J.; Cheeseman, J. R.; Frisch, M. J. J. Phys. Chem. A 2001, 105, 5356–5371. Kaminsky´, J.; Tomcˇa´kova´, K.; Bourˇ, P.; Raich, I. J. Comput. Chem. 2010, 31, 2213–2224. Sadlej, A. J. Collect. Czech. Chem. Commun. 1988, 53, 1995–2016. Ruud, K.; Stephens, P. J.; Devlin, F. J.; et al. Chem. Phys. Lett. 2003, 373, 606–614. Campos, C. T.; Jorge, F. E.; Silva, T. P.; Coppo, M. R. Chem. Phys. Lett. 2010, 494, 170–173. Autschbach, J.; Ziegler, T.; van Gisbergen, S. J. A.; Baerends, E. J. J. Chem. Phys. 2002, 116, 6930–6940. Ruud, K.; Helgaker, T. Chem. Phys. Lett. 2002, 352, 533–539. Kundrat, M. D.; Autschbach, J. J. Phys. Chem. A 2006, 110, 4115–4123. Autschbach, J. Chirality 2009, 21, E116–E152.

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

543

34. Cappelli, C.; Corni, S.; Mennucci, B.; Cammi, R.; Tomasi, J. J. Phys. Chem. A 2002, 106, 12331–12339; Aamouche, A.; Devlin, F. J.; Stephens, P. J. J. Am. Chem. Soc. 2000, 122, 7358–7367. 35. Kundrat, M. D.; Autschbach, J. J. Chem. Theory Comput. 2008, 4, 1902–1914. 36. Mukhopadhyay, P.; Zuber, G.; Goldsmith, M. R.; Wipf, P.; Beratan, D. N. Chem. Phys. Lett. 2006, 7, 2483–2486; Dracˇ´ınsky´, M.; Kaminsky´, J.; Bourˇ, P. J. Phys. Chem. B 2009, 113, 14698–14707. 37. Goldsmith, M. R.; Jayasuria, N.; Beratan, D. N.; Wipf, P. J. Am. Chem. Soc. 2003, 125, 15696–15697. 38. Mukhopadhyay, P.; Wipf, P.; Beratan, D. N. Acc. Chem. Res. 2009, 42, 809–819. 39. Pedersen, T. B.; Kongsted, J.; Crawford, T. D.; Ruud, K. J. Chem. Phys. 2009, 130, 034310; Ruud, K.; Zanasi, R. Angew. Chem. Int. Ed. 2005, 44, 3594–3596. 40. Kongsted, J.; Ruud, K. Chem. Phys. Lett. 2008, 451, 226–232. 41. Wiberg, K. B.; Vaccaro, P. H.; Cheeseman, J. R. J. Am. Chem. Soc. 2003, 125, 1888–1896; Wilson, S. M.; Wiberg, K. B.; Cheeseman, J. R.; Frisch, M. J.; Vaccaro, P. H. J. Phys. Chem. A 2005, 109, 11752–11764. 42. Barron, L. D. Molecular Light Scattering and Optical Activity; Cambridge University Press: Cambridge, 2004. 43. Moscowitz, A. Tetrahedron 1961, 13, 48–56. 44. Polavarapu, P. L. J. Phys. Chem. A 2005, 109, 7013–7023; Polavarapu, P. L. Chirality 2006, 18, 348–356. 45. DeVoe, H. J. Chem. Phys. 1965, 43, 3199–3208. 46. Cech, C. L.; Hug, W.; Tinoco, I. J. Biopolymers 1976, 15, 131–152. 47. Bayley, P. M.; Nielsen, E. B.; Schellman, J. A. J. Phys. Chem. 1969, 73, 228–243. 48. Bulheller, B. M.; Rodger, A.; Hirst, J. D. Phys. Chem. Chem. Phys. 2007, 9, 2020–2035. 49. Bulheller, B. M.; Hirst, J. D. Bioinformatics 2009, 25, 539–540. 50. Foresman, J. B.; Head-Gordon, M.; Pople, J. A.; Frisch, M. J. J. Phys. Chem. 1992, 96, 135–149. 51. Olsen, J.; Roos, B. O.; Jorgensen, P.; Jensen, H. J. J. Chem. Phys. 1988, 89, 2185–2192. 52. Andersson, K.; Malmqvist, P. A˚.; Roos, B. O. J. Chem. Phys. 1992, 96, 1218–1226. 53. Schleyer, P. R.; Allinger, N. L.; Clark, T.; et al. The Encyclopedia of Computational Chemistry; Wiley: Chichester, UK, 1998. 54. Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029–7039. 55. Furche, F.; Alhrichs, R. J. Chem. Phys. 2004, 121, 12772–12773. 56. Kotzian, M.; Rosch, N.; Zerner, M. C. Theor. Chim. Acta 1992, 81, 201–222. 57. Roos, B. O. Int. J. Quantum Chem. 1980, S14, 175–189. 58. Cˇ´ızˇek, J. On the Use of the Cluster Expansion and the Technique of Diagrams in Calculations of Correlation Effects in Atoms and Molecules. In Advances in Chemical Physics: Correlation Effects in Atoms and Molecules; Hariharan, P. C., Ed.; Wiley Interscience: New York, 1969, pp 35–91, Bartlett, R. J. Modern Ideas in Coupled Cluster Methods; World Scientific: Singapore, 1997. 59. Das, A.; Hasegawa, J.; Miyahara, T.; Ehara, M.; Nakatsuji, H. J. Comput. Chem. 2003, 24, 1421–1431. 60. Christiansen, O.; Koch, H.; Jørgensen, P. J. Chem. Phys. 1995, 103, 7429. 61. Furche, F.; Burke, K. Ann. Rep. Comput. Chem. 2005, 1, 19–30; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2002, 117, 7433–7447. 62. Autschbach, J.; Ziegler, T. Coord. Chem. Rev. 2003, 238–239, 83–126. 63. Santoro, F.; Barone, V. Int. J. Quantum Chem. 2009, 110, 476–486. 64. Furche, F.; Ahlrichs, R.; Wachsmann, C.; et al. J. Am. Chem. Soc. 2000, 122, 1717–1724. 65. Julı´nek, O.; Setnicˇka, V.; Mikla´sˇova´, N.; et al. J. Phys. Chem. A 2009, 113, 10717–10725. 66. Sreerama, N.; Manning, M. C.; Powers, M. E.; et al. Biochemistry 1999, 38, 10814–10822. 67. Pancˇosˇka, P.; Bitto, E.; Janota, V.; et al. Protein Sci. 1995, 4, 1384–1401. 68. Kelly, S. M.; Jess, T. J.; Price, N. C. Biochim. Biophys. Acta 2005, 1751, 119–139. 69. Kaminsky´, J.; Kubelka, J.; Bourˇ, P. J. Phys. Chem. A 2011, 115, 1734–1742. 70. Sˇebek, J.; Bourˇ, P. J. Phys. Chem. A 2008, 112, 2920–2929. 71. Julı´nek, O.; Krupicˇka, M.; Lindner, W.; Urbanova´, M. Phys. Chem. Chem. Phys. 2010, 12, 11487–11497. 72. Keiderling, T. A. J. Chem. Phys. 1981, 75, 3639–3641. 73. Funk, T.; Deb, A.; George, S. J.; Wang, H. X.; Cramer, S. P. Coord. Chem. Rev. 2005, 249, 3–30. 74. Stephens, P. J. J. Chem. Phys. 1970, 52, 3489–3516. 75. Seth, M.; Ziegler, T.; Autschbach, J. The J. Chem. Phys. 2008, 129, 104105. 76. Michl, J. J. Am. Chem. Soc. 1978, 100, 6801–6811. 77. Solheim, H.; Ruud, K.; Coriani, S.; Norman, P. J. Chem. Phys. 2008, 128, 094103. 78. Lowry, T. M. Optical Rotatory Power; Longmans, Green and Co.: London, 1935. 79. Hsu, E. C.; Holzwarth, G. J. Chem. Phys. 1973, 59, 4678–4685. 80. Holzwarth, G.; Hsu, E. C.; Mosher, H. S.; Faulkner, T. R.; Moscowitz, A. J. Am. Chem. Soc. 1974, 96, 251–252. 81. Nafie, L. A.; Keiderling, T. A.; Stephens, P. J. J. Am. Chem. Soc. 1976, 98, 2715–2723. 82. Polyanichko, A. M.; Andrushchenko, V.; Bourˇ, P.; Wieser, H. Vibrational Circular Dichroism Studies of Biological Macromolecules and their Complexes. In Circular Dichroism: Theory and Spectroscopy; Rodgers, D. S., Ed.; Nova: Hauppauge, NY, 2011. 83. Schellman, J. A. J. Chem. Phys. 1973, 58, 2882–2886. 84. Faulkner, T. R.; Moscowitz, A.; Holzwarth, G.; Hsu, E. C.; Mosher, H. S. J. Am. Chem. Soc. 1974, 96, 252–253. 85. Holzwarth, G.; Chabay, I. J. Chem. Phys. 1972, 57, 1632–1635. 86. Bourˇ, P.; Keiderling, T. A. J. Am. Chem. Soc. 1992, 114, 9100–9105. 87. Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. 1990, 141, 189–196. 88. Coriani, S.; Thorvaldsen, A. J.; Kristensen, K.; Jørgensen, P. Phys. Chem. Chem. Phys. 2011, 13, 4224–4229. 89. Yang, D.; Rauk, A. J. J. Chem. Phys. 1992, 97, 6517–6534. 90. Bourˇ, P.; McCann, J.; Wieser, H. J. Chem. Phys. 1998, 108, 8782–8789. 91. Stephens, P. J. J. Phys. Chem. 1987, 91, 1712–1715. 92. Cheeseman, J. R.; Frisch, M. J.; Devlin, F. J.; Stephens, P. J. Chem. Phys. Lett. 1996, 252, 211–220. 93. Bak, K. L.; Bludsky´, O.; Jorgensen, P. J. Chem. Phys. 1995, 103, 10548–10555. 94. Urbanova´, M.; Setnicˇka, V.; Devlin, F.; Stephens, P. J. J. Am. Chem. Soc. 2005, 127, 6700–6711. 95. Daneˇcˇek, P.; Kapita´n, J.; Baumruk, V.; et al. J. Chem. Phys. 2007, 126, 224513. 96. Bourˇ, P.; Keiderling, T. A. J. Am. Chem. Soc. 1993, 115, 9602–9607. 97. Bourˇ, P.; Sopkova´, J.; Bedna´rova´, L.; Malonˇ, P.; Keiderling, T. A. J. Comput. Chem. 1997, 18, 646–659. 98. Bourˇ, P.; Keiderling, T. A. J. Chem. Phys. 2002, 117, 4126–4132. 99. Andrushchenko, V.; Wieser, H.; Bourˇ, P. J. Phys. Chem. B 2002, 106, 12623–12634. 100. Bourˇ, P.; Michalı´k, D.; Kapita´n, J. J. Chem. Phys. 2005, 122, 144501.

544

101. 102. 103. 104. 105. 106. 107. 108. 109. 110. 111. 112. 113. 114. 115. 116. 117. 118. 119. 120. 121. 122.

Spectroscopic Analysis: Ab initio Calculation of Chiroptical Spectra

Kwac, K.; Lee, K. K.; Han, J.; Oh, K. I.; Cho, M. J. Chem. Phys. 2008, 128, 105106. Devine, T. R.; Keiderling, T. A. J. Phys. Chem. 1984, 88, 390–394. Bourˇ, P.; Tam, C. N.; Wang, B.; Keiderling, T. A. Mol. Phys. 1996, 87, 299–318. Barron, L. D.; Bogaard, M. P.; Buckingham, A. D. J. Am. Chem. Soc. 1973, 95, 603–605. Nafie, L. A.; Che, D. Theory and Measurement of Raman Optical Activity. In Modern Nonlinear Optics, Part 3; Evans, M., Kielich, S., Eds.; Wiley: New York, 1994, Vol. 85; pp 105–149. Ruud, K.; Thorvaldsen, J. Chirality 2009, 21, E54–E67. Angeli, C.; Bak, K. L.; Bakken, V.; et al. Dalton, A Molecular Electronic Structure Program, Release 2.0; University of Oslo: Oslo, 2009. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. Gaussian 09, Revision A.02; Gaussian: Wallingford, CT, 2009. Ahlrichs, R.; Bar, M.; Baron, H.-P.; et al. Turbomole, Version 5; Quantum Chemistry Group, University of Karlsruhe: Karlsruhe, 1998. Bourˇ, P.; Baumruk, V.; Hanzlı´kova´, J. Collect. Czech. Chem. Commun. 1997, 62, 1384–1395. Kapita´n, J.; Baumruk, V.; Kopecky´, V., Jr.; Pohl, R.; Bourˇ, P. J. Am. Chem. Soc. 2006, 128, 13451–13462. Helgaker, T.; Ruud, K.; Bak, K. L.; Joergensen, P.; Olsen, J. Faraday Discuss. 1994, 99, 165–180. Choi, J. H.; Cho, M. H. J. Chem. Phys. 2009, 130, 014503. Luber, S.; Reiher, M. Chem. Phys. Chem. 2009, 10, 2049–2057. Polavarapu, P. L. Angew. Chem. Int. Ed. 2002, 41, 4544–4546. Haesler, J.; Schindelholz, I.; Riguet, E.; Bochet, C. G.; Hug, W. Nature 2007, 446, 526–529. Kapita´n, J.; Baumruk, V.; Kopecky´, V., Jr.; Bourˇ, P. J. Phys. Chem. A 2006, 110, 4689–4696. Luber, S.; Reiher, M. J. Phys. Chem. B 2010, 114, 1057–1063. Yamamoto, S.; Straka, M.; Watarai, H.; Bourˇ, P. Phys. Chem. Chem. Phys. 2010, 12, 11021–11032. Budeˇsˇ´ınsky´, M.; Daneˇcˇek, P.; Bedna´rova´, L.; et al. J. Phys. Chem. A 2008, 112, 8633–8640. ADF 2007. SCM, Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands, 2007, http://www.scm.com. Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2011, 7, 291–309.