Classical Magnetic Dipole Moments for the ... - ACS Publications

Jan 15, 2016 - spectra are given by the Fourier transform of the cross-correlation function of magnetic dipole moment and electric dipole moment. We o...
0 downloads 0 Views 666KB Size
Download PDF
Letter pubs.acs.org/JPCL

Classical Magnetic Dipole Moments for the Simulation of Vibrational Circular Dichroism by ab Initio Molecular Dynamics Martin Thomas* and Barbara Kirchner* Mulliken Center for Theoretical Chemistry, Rheinische Friedrich-Wilhelms-Universität Bonn, Beringstraße 4, 53115 Bonn, Germany S Supporting Information *

ABSTRACT: We present a new approach for calculating vibrational circular dichroism spectra by ab initio molecular dynamics. In the context of molecular dynamics, these spectra are given by the Fourier transform of the cross-correlation function of magnetic dipole moment and electric dipole moment. We obtain the magnetic dipole moment from the electric current density according to the classical definition. The electric current density is computed by solving a partial differential equation derived from the continuity equation and the condition that eddy currents should be absent. In combination with a radical Voronoi tessellation, this yields an individual magnetic dipole moment for each molecule in a bulk phase simulation. Using the chiral alcohol 2-butanol as an example, we show that experimental spectra are reproduced very well. Our approach requires knowing only the electron density in each simulation step, and it is not restricted to any particular electronic structure method.

V

is the computation of the magnetic dipole moment. In ref 20, the MD force field was parametrized with atomic axial tensors from static quantum chemistry in terms of magnetic field perturbation theory. The authors of refs 11 and 12 are directly aiming for a combination of nuclear velocity perturbation theory with ab initio molecular dynamics (AIMD). Another path was taken in refs 19, 21, and 22, in which the definition of the classical magnetic dipole moment caused by point charges was used. The point charges on the atoms were either set to fixed values,19 or they were obtained by population analyses in a quantum mechanical/molecular mechanical setup. 21,22 Although a good agreement with the experiment was achieved, it should be noted that atomic partial charges cannot be defined in a unique way, and the result always depends on the particular population analysis that is selected.22 In this work, we present a novel approach for calculating the magnetic dipole moment in an AIMD simulation. We follow the previous idea to employ the classical definition of the magnetic dipole moment, but we do not rely on atomic partial charges. Instead, we process the complete information about the electron density that is inherently available in AIMD. In combination with our recently developed technique of applying a radical Voronoi tessellation to the electron density,23 this allows a straightforward calculation of an individual magnetic dipole moment for each molecule in a bulk phase system. After the explanation of our method to calculate magnetic dipole moments in the following paragraphs, we show the application to the small chiral alcohol 2-butanol as an example system. We investigate the differences between gas phase and bulk phase,

ibrational circular dichroism (VCD) is an important analytical technique for studying the vibrational spectra of chiral molecules.1−3 It measures the difference in the absorption of left and right circularly polarized infrared (IR) radiation. The two enantiomers of a chiral molecule show the same VCD intensities but with opposite sign, so their spectra are mirror images with respect to the wavenumber axis. Thus, VCD can be used to determine the absolute configuration of a chiral molecule. The ab initio prediction of VCD spectra is an interesting subject of theoretical chemistry. Some quantum chemistry software packages provide the possibility of obtaining VCD intensities by static calculations on the basis of density functional theory. These implementations rely either on magnetic field perturbation theory4−7 or on nuclear velocity perturbation theory8−12 to compute atomic axial tensors and rotational strengths. To reproduce experimental conditions, it is often necessary to include the influence of the molecular surrounding in the bulk phase or in a solution. In the context of static quantum chemistry, this is possible by implicit solvation models and microsolvation approaches with explicit solvent molecules, as has been done, for example, to study the VCD spectra of several biomolecules.13−16 Such solute−solvent clusters can be extracted from molecular dynamics (MD) simulations,17,18 but it would be desirable to directly employ MD without modifications of the molecular surrounding to capture the dynamic nature of the intermolecular interaction network and to cover the conformational flexibility in the bulk phase. It has been shown that the VCD spectrum can be obtained from an MD simulation by Fourier transforming the crosscorrelation function of the electric dipole moment and the magnetic dipole moment.19−21 The key point in this approach © XXXX American Chemical Society

Received: December 11, 2015 Accepted: January 15, 2016

509

DOI: 10.1021/acs.jpclett.5b02752 J. Phys. Chem. Lett. 2016, 7, 509−513

Letter

The Journal of Physical Chemistry Letters and the simulated spectra are compared to experimental data. Simulations of propylene oxide and α-pinene as further examples can be found in the Supporting Information. Based on an MD simulation, the VCD spectrum can be calculated by Fourier transforming the time cross-correlation function of the electric dipole moment, μ̇ (t), and the magnetic dipole moment, m(t).19−21 Following our previous expression for IR spectra,23,24 the VCD spectrum is defined as ΔA(ν)̃ =

This ansatz ensures that the electric current follows an irrotational vector field and flows only in spatial regions with a finite electron density. The same model is used, for example, for a potential flow in fluid dynamics, which is employed to describe the irrotational flow of a homogeneous fluid with negligible frictional forces. This is a reasonable analogon of the molecular electron density in the classical limit. Assuming the absence of eddy currents is certainly an approximation, but it works very well for the systems selected in this work. Nevertheless, a comparison to a proper quantum calculation of the electric current density could be subject of future work. Inserting ansatz 4 into the continuity equation (eq 3) leads to the partial differential equation

∞

NA

∫ [⟨μ̇ (τ)·ṁ (τ + t )⟩τ 12ε c k T −∞ 2

0

B

− ⟨ṁ (τ ) ·μ̇ (τ + t )⟩τ ]exp( −2π icνt̃ )dt

(1)

where ν̃ denotes the wavenumber and t and τ denote the time; NA is the Avogadro constant, ε0 the vacuum permittivity, c the speed of light in vacuum, kB the Boltzmann constant, and T the average simulation temperature; and ⟨μ̇ (τ)·ṁ (τ+ t)⟩τ denotes the cross-correlation function of the time derivatives of the electric dipole moment and the magnetic dipole moment. The key point in this approach is the calculation of the magnetic dipole moment, for which we propose a novel method as described in the following. We have recently shown that molecular electric dipole moments in a bulk phase AIMD simulation can be computed by a radical Voronoi tessellation and subsequent integration of the electron density according to the classical definition.23 The corresponding classical expression for the magnetic dipole moment reads as m=

1 2

∫ r × j(r)dr

∂ρ(r, t ) = ∇ρ(r, t ) ·∇α(r, t ) + ρ(r, t )Δα(r, t ) ∂t

which has to be solved to find the scalar field α(r, t) that finally leads to the electric current density j(r, t). For a numerical solution, eq 5 is transformed to a system of linear equations by a finite difference approach. The discretization grid is naturally given, for example, by the grid of the Gaussian cube file of the electron density data written by the simulation software. This leads to a nonsymmetric equation system with a size in the order of 106 equations and unknowns for the systems studied in this work. It can be solved iteratively, for example, by the BiCGstab(l) algorithm,26−28 which we have implemented in TRAVIS.23,24,29,30 Solving eq 5 requires a certain computational effort, but in combination with AIMD on the basis of density functional theory, it constitutes only a minor contribution to the total resources needed to simulate a VCD spectrum (see the Supporting Information for further details). The definition of the magnetic dipole moment according to eq 2 depends on the choice of the coordinate system. This is the same effect that occurs for the electric dipole moment of a charged system. We have discussed before23 that the most reasonable choice for the reference point in the context of vibrational spectra is the molecular center of mass because it remains fixed in all normal modes; therefore, it is independent of the internal degrees of freedom. For the same reason, we use the center of mass as reference point for molecular magnetic dipole moments here (see the Supporting Information for further details). This choice of a different reference point for each molecule under periodic boundary conditions is another approximation, but it works very well for the systems studied in this work. To show the applicability of our approach, we carried out an AIMD simulation of (R)-2-butanol in the gas phase using density functional theory with the BLYP exchange-correlation functional31,32 and a double-ζ basis set (see the Supporting Information for further details). Because not only the magnetic dipole moment but also the electric dipole moment enters the calculation of the VCD spectrum, very little effort is necessary to obtain the IR spectrum too. In Figure 1, both spectra are compared to experimental data of a dilute solution of (R)-2butanol in carbon disulfide,33 which should provide a reasonable approximation to the gas phase.34 Due to the restrictions imposed by the experimental setup and the IR absorption of the solvent, the measured spectra are available only between 900 and 1400 cm−1, but this is the most important wavenumber region for 2-butanol. Apart from the slight redshift of many bands, which we have observed in other systems when applying the same electronic

(2)

with the electric current density j(r) and the spatial coordinate r. The nuclear contribution to the magnetic dipole moment can straightforwardly be obtained according to the definition for point charges as mentioned above, replacing the atomic partial charges by the nuclear core charges. As soon as the electric current density caused by the electrons is known, individual molecular magnetic dipole moments can be calculated in the bulk phase completely analogous to electric dipole moments, leading to a straightforward calculation of individual molecular VCD spectra in such systems. It should be noted that other partitioning schemes for the electron density such as, for example, Bader’s atoms in molecules25 could be used for the same purpose. To obtain the electric current density caused by the electrons, we start with the continuity equation ∂ρ(r, t ) + ∇·j(r, t ) = 0 ∂t

(3)

which means that changes in the electron density, ρ(r, t), act as sources and sinks of an electric current. The time derivative can be approximated by saving the electron density in each simulation step and calculating finite differences afterward. The continuity equation does not uniquely determine the electric current density because the divergence is not changed by adding an arbitrary solenoidal vector field. In other words, arbitrary eddy currents do not appear as changes of the electron density. Thus, we introduce as further condition that the electric current density should be the product of the electron density and a conservative velocity field, v(r, t) = −∇α(r, t): j(r, t ) = −ρ(r, t )∇α(r, t )

(5)

(4) 510

DOI: 10.1021/acs.jpclett.5b02752 J. Phys. Chem. Lett. 2016, 7, 509−513

Letter

The Journal of Physical Chemistry Letters

Figure 1. Simulated IR and VCD spectra of (R)-2-butanol in the gas phase. Experimental IR and VCD spectra of a 0.029 M (R)-2-butanol solution in carbon disulfide are reproduced from ref 33. Copyright 2000 American Chemical Society.

Figure 2. Simulated IR and VCD spectra of (R)-2-butanol and (S)-2butanol in the bulk phase. The experimental IR and VCD spectra of the (R)-enantiomer are reproduced from ref 33. Copyright 2000 American Chemical Society.

structure method before,23,24,30 the IR and VCD spectra are in good agreement. In particular, the experimentally observed VCD pattern with a positive peak at 1076 cm−1 followed by a strong negative peak at 1143 cm−1 and another positive peak at 1241 cm−1 appears in the same way in the simulated spectrum. Although the positions of the bands are slightly shifted to 1078, 1118, and 1227 cm−1 in the AIMD, their intensity ratios are reproduced very well. The coincidence of simulation and experiment is of particular interest regarding the conformational flexibility of 2-butanol. Both the C−O bond and the central C−C bond have three distinct torsional minima on the potential energy surface, resulting in nine conformations of the molecule. Static calculations reveal that these possess significantly different VCD spectra on their own,33,34 but it is apparent that the AIMD simulation properly averages over all important conformers. Due to the ability to form intermolecular hydrogen bonds, significant differences in the spectra of the gas phase and the liquid phase of 2-butanol are expected. Therefore, we performed AIMD simulations of 16 2-butanol molecules under periodic boundary conditions (see the Supporting Information for details). We carried out two independent simulations, one for the (R)-configuration and one for the (S)configuration, to check if our approach results in mirror symmetric VCD spectra for the two enantiomers. The obtained IR and VCD spectra are compared to experimental data33 in Figure 2. The IR spectrum is shown only for the (R)enantiomer because the one of the (S)-enantiomer is almost equal, as expected. Minor differences occur because of the limited system size and the finite length of the trajectories. Apart from the aforementioned redshift of many fingerprint

bands, the AIMD simulation reproduces the effect of the hydrogen bonding on the IR spectrum very well, because the broadening of the OH stretching band above 3000 cm−1 and the appearance of the librations as another broad band around 600 cm−1 are in good agreement with the experiment. The simulated VCD spectra show the expected result, because each band with a positive sign in the (R)-configuration has a negative sign in the (S)-configuration and vice versa. Thus, our model is valid regarding this basic property of VCD spectroscopy. Minor differences in the intensity ratios are again caused by the limited system size and the finite length of the trajectories. Furthermore, the pattern of positive and negative peaks in the experimental VCD spectrum of the (R)enantiomer is reproduced very well by the AIMD. The only exception is the positive band found at 968 cm−1 in the experiment, which is strongly underestimated in the simulation. However, also the IR spectrum is not fully consistent with the experiment in this wavenumber region, so this is unlikely to be a particular issue of our VCD approach. When gas phase and liquid phase are compared, the most important change is the vanishing positive peak at 1241 cm−1. According to a normal coordinate analysis of the AIMD trajectories,30 this band is connected to the COH bending mode, which is shifted to the region of the CH bending modes around 1400 cm−1 in the liquid. The behavior of this vibration due to the hydrogen bonding in the bulk phase is reproduced very well by the AIMD simulations. In summary, we presented a novel method of calculating molecular magnetic dipole moments and VCD spectra by AIMD simulations. The electron density inherently available in AIMD is processed and a partial differential equation is solved 511

DOI: 10.1021/acs.jpclett.5b02752 J. Phys. Chem. Lett. 2016, 7, 509−513

Letter

The Journal of Physical Chemistry Letters

(4) Stephens, P. J. Theory of Vibrational Circular Dichroism. J. Phys. Chem. 1985, 89, 748−752. (5) Stephens, P. J. Gauge Dependence of Vibrational Magnetic Dipole Transition Moments and Rotational Strengths. J. Phys. Chem. 1987, 91, 1712−1715. (6) Cheeseman, J.; Frisch, M.; Devlin, F.; Stephens, P. Ab Initio Calculation of Atomic Axial Tensors and Vibrational Rotational Strengths Using Density Functional Theory. Chem. Phys. Lett. 1996, 252, 211−220. (7) Nicu, V. P.; Neugebauer, J.; Wolff, S. K.; Baerends, E. J. A vIbrational Circular Dichroism Implementation within a Slater-TypeOrbital Based Density Functional Framework and its Application to Hexa- and Hepta-Helicenes. Theor. Chem. Acc. 2008, 119, 245−263. (8) Nafie, L. A.; Freedman, T. B. Vibronic Coupling Theory of Infrared Vibrational Transitions. J. Chem. Phys. 1983, 78, 7108−7116. (9) Nafie, L. A. Adiabatic Molecular Properties beyond the BornOppenheimer Approximation. Complete Adiabatic Wave Functions and Vibrationally Induced Electronic Current Density. J. Chem. Phys. 1983, 79, 4950−4957. (10) Buckingham, A.; Fowler, P.; Galwas, P. Velocity-Dependent Property Surfaces and the Theory of Vibrational Circular Dichroism. Chem. Phys. 1987, 112, 1−14. (11) Scherrer, A.; Vuilleumier, R.; Sebastiani, D. Nuclear Velocity Perturbation Theory of Vibrational Circular Dichroism. J. Chem. Theory Comput. 2013, 9, 5305−5312. (12) Scherrer, A.; Agostini, F.; Sebastiani, D.; Gross, E. K. U.; Vuilleumier, R. Nuclear Velocity Perturbation Theory for Vibrational Circular Dichroism: An Approach Based on the Exact Factorization of the Electron-Nuclear Wave Function. J. Chem. Phys. 2015, 143, 074106. (13) Poopari, M. R.; Dezhahang, Z.; Yang, G.; Xu, Y. Conformational Distributions of N-Acetyl-L-cysteine in Aqueous Solutions: A Combined Implicit and Explicit Solvation Treatment of VA and VCD Spectra. ChemPhysChem 2012, 13, 2310−2321. (14) Cappelli, C.; Bloino, J.; Lipparini, F.; Barone, V. Toward Ab Initio Anharmonic Vibrational Circular Dichroism Spectra in the Condensed Phase. J. Phys. Chem. Lett. 2012, 3, 1766−1773. (15) Rode, J. E.; Dobrowolski, J. C.; Sadlej, J. Prediction of LMethionine VCD Spectra in the Gas Phase and Water Solution. J. Phys. Chem. B 2013, 117, 14202−14214. (16) Orestes, E.; Bistafa, C.; Rivelino, R.; Canuto, S. Including Thermal Disorder of Hydrogen Bonding to Describe the Vibrational Circular Dichroism Spectrum of Zwitterionic L-Alanine in Water. J. Phys. Chem. A 2015, 119, 5099−5106. (17) Kundrat, M. D.; Autschbach, J. Ab Initio and Density Functional Theory Modeling of the Chiroptical Response of Glycine and Alanine in Solution Using Explicit Solvation and Molecular Dynamics. J. Chem. Theory Comput. 2008, 4, 1902−1914. (18) Tedesco, D.; Zanasi, R.; Kirchner, B.; Bertucci, C. Short-Range Solvation Effects on Chiroptical Properties: A Time-Dependent Density Functional Theory and Ab Initio Molecular Dynamics Computational Case Study on Austdiol. J. Phys. Chem. A 2014, 118, 11751−11757. (19) Abbate, S.; Longhi, G.; Kwon, K.; Moscowitz, A. The Use of Cross-Correlation Functions in the Analysis of Circular Dichroism Spectra. J. Chem. Phys. 1998, 108, 50−62. (20) Horníček, J.; Kaprálová, P.; Bouř, P. Simulations of Vibrational Spectra from Classical Trajectories: Calibration with Ab Initio Force Fields. J. Chem. Phys. 2007, 127, 084502. (21) Yang, S.; Cho, M. Direct Calculations of Vibrational Absorption and Circular Dichroism Spectra of Alanine Dipeptide Analog in Water: Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations. J. Chem. Phys. 2009, 131, 135102. (22) Choi, J.-H.; Cho, M. Direct Calculations of Mid- and Near-IR Absorption and Circular Dichroism Spectra of Chiral Molecules Using QM/MM Molecular Dynamics Simulation Method. J. Chem. Theory Comput. 2011, 7, 4097−4103.

to obtain the electric current density on the basis of the continuity equation and the condition that there should not be any eddy currents. According to the classical definition and in combination with a radical Voronoi tessellation, this provides an individual magnetic dipole moment for each molecule in the simulation cell. The VCD spectrum is given by the Fourier transform of the cross-correlation function of magnetic dipole moment and electric dipole moment. For the example system 2-butanol, our approach yields VCD spectra in very good agreement with experimental data. All essential features of the spectra are reproduced by the simulations, in particular, also differences between the gas phase and the bulk phase are modeled very well. Two further examples, propylene oxide and α-pinene, which confirm the conclusions can be found in the Supporting Information. Our proposed method of computing magnetic dipole moments directly relies on the electron density from the AIMD simulation. It is not required to use a population analysis to obtain atomic partial charges, and furthermore, no modifications of the electronic structure code are necessary. Thus, the approach can readily be combined with any electronic structure method that is able to provide the electron density and it is not restricted to density functional theory with the GGA functional BLYP31,32 and a double-ζ basis set as applied in this work. This turns it into a valuable general tool for the simulation of VCD spectra.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.5b02752. Computational details, an investigation of the convergence of the magnetic dipole moments with respect to several numerical parameters that enter their calculation, a study of the influence of the reference point for the magnetic dipole moments, and the results for propylene oxide and α-pinene as further example systems (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*Phone: +49 228 7360439. Fax: +49 228 739064. E-mail: [email protected]. *Phone: +49 228 7360439. Fax: +49 228 739064. E-mail: [email protected] Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Stefan Grimme and Christoph Jacob for helpful discussions. Financial support by the DFG through the SFB 813 is gratefully acknowledged.

■

REFERENCES

(1) Nafie, L. A. Infrared and Raman Vibrational Optical Activity: Theoretical and Experimental Aspects. Annu. Rev. Phys. Chem. 1997, 48, 357−386. (2) Magyarfalvi, G.; Tarczay, G.; Vass, E. Vibrational Circular Dichroism. WIREs Comput. Mol. Sci. 2011, 1, 403−425. (3) Wu, T.; You, X.-Z.; Bouř, P. Applications of Chiroptical Spectroscopy to Coordination Compounds. Coord. Chem. Rev. 2015, 284, 1−18. 512

DOI: 10.1021/acs.jpclett.5b02752 J. Phys. Chem. Lett. 2016, 7, 509−513

Letter

The Journal of Physical Chemistry Letters (23) Thomas, M.; Brehm, M.; Kirchner, B. Voronoi Dipole Moments for the Simulation of Bulk Phase Vibrational Spectra. Phys. Chem. Chem. Phys. 2015, 17, 3207−3213. (24) Thomas, M.; Brehm, M.; Fligg, R.; Vöhringer, P.; Kirchner, B. Computing Vibrational Spectra from Ab Initio Molecular Dynamics. Phys. Chem. Chem. Phys. 2013, 15, 6608−6622. (25) Bader, R. F. W. A Quantum Theory of Molecular Structure and its Applications. Chem. Rev. 1991, 91, 893−928. (26) van der Vorst, H. Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM J. Sci. Stat. Comput. 1992, 13, 631−644. (27) Sleijpen, G. L.; Fokkema, D. R. BiCGstab(l) for Linear Equations Involving Unsymmetric Matrices with Complex Spectrum. Electronic Transactions on Numerical Analysis 1993, 1, 11−32. (28) Sleijpen, G.; van der Vorst, H. Reliable Updated Residuals in Hybrid Bi-CG Methods. Computing 1996, 56, 141−163. (29) Brehm, M.; Kirchner, B. TRAVIS - A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories. J. Chem. Inf. Model. 2011, 51, 2007−2023. (30) Thomas, M.; Brehm, M.; Hollóczki, O.; Kelemen, Z.; Nyulászi, L.; Pasinszki, T.; Kirchner, B. Simulating the Vibrational Spectra of Ionic Liquid Systems: 1-Ethyl-3-methylimidazolium Acetate and its Mixtures. J. Chem. Phys. 2014, 141, 024510. (31) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (32) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (33) Wang, F.; Polavarapu, P. L. Vibrational Circular Dichroism: Predominant Conformations and Intermolecular Interactions in (R)(−)-2-Butanol. J. Phys. Chem. A 2000, 104, 10683−10687. (34) Shin, S.; Nakata, M.; Hamada, Y. Analysis of Vibrational Circular Dichroism Spectra of (S)-(+)-2-Butanol by Rotational Strengths Expressed in Local Symmetry Coordinates. J. Phys. Chem. A 2006, 110, 2122−2129.

513

DOI: 10.1021/acs.jpclett.5b02752 J. Phys. Chem. Lett. 2016, 7, 509−513