Direct MD Simulations of Terahertz Absorption and 2D Spectroscopy

Feb 10, 2014 - A direct molecular dynamics simulation of the THz spectrum of a molecular crystal is presented. A time-dependent electric field is adde...
0 downloads 0 Views 2MB Size
Download PDF
Letter pubs.acs.org/JPCL

Direct MD Simulations of Terahertz Absorption and 2D Spectroscopy Applied to Explosive Crystals G. Katz,† S. Zybin,‡ W. A. Goddard, III,‡ Y. Zeiri,¶,§ and R. Kosloff*,† †

Fritz Haber Research Center for Molecular Dynamics, Hebrew University of Jerusalem, Jerusalem 91904, Israel Materials and Process Simulation Center, California Institute of Technology, Pasadena, California 91125, United States ¶ Chemistry Department, NRCN, P.O. Box 9001, Beer-Sheva 84190, Israel § Bio-Medical Engineering, Ben-Gurion University, Beer-Sheva 84105, Israel ‡

S Supporting Information *

ABSTRACT: A direct molecular dynamics simulation of the THz spectrum of a molecular crystal is presented. A time-dependent electric field is added to a molecular dynamics simulation of a crystal slab. The absorption spectrum is composed from the energy dissipated calculated from a series of applied pulses characterized by a carrier frequency. The spectrum of crystalline cyclotrimethylenetrinitramine (RDX) and triacetone triperoxide (TATP) were simulated with the ReaxFF force field. The proposed direct method avoids the linear response and harmonic approximations. A multidimensional extension of the spectroscopy is suggested and simulated based on the nonlinear response to a single polarized pulse of radiation in the perpendicular polarization direction. SECTION: Spectroscopy, Photochemistry, and Excited States

R

spectrum that hinders the task of assignment. Can computational tools supply a reliable simulation that can be used for the interpretation task in the THz band? For the higher vibrational frequencies, quantum chemistry codes have become essential. The first step is geometry optimization, which is subsequently followed by normal-mode analysis. Combined with the ability of quantum chemistry to calculate dipole functions, an assignment is suggested22 supplemented with the addition of anharmonic corrections.23 The extension of these procedures to the THz regime is difficult.14,24−26 Typically, one is dealing with very large molecules or extended molecular systems. The soft modes considered are strongly coupled and therefore show strong anharmonic features. To overcome the problem of simulating the spectrum, molecular dynamics (MD) techniques have been employed to assign normal modes.15,27,28 Neglecting quantum effects is justified because at the THz spectral range in the temperatures considered, tunnelling splitting and zero-point motion are suppressed. The present study is aimed at the development of a computational tool for simulating and assigning THz spectroscopy directly. Such a method should not rely on linear response and the harmonic approximation. For the simulation, we employ a MD approach applied to a crystal slab with periodic boundary conditions. The force field used is ReaxFF, which is

emote sensing is based on the ability of light to interact with the subject matter and convey the information to the detector. Different light−matter interactions can be selected by a choice of frequency band.1,2 The terahertz band (THz) has become a candidate for remote sensing due to its ability to penetrate dry material3−6 and the development of new THz sources.7 The remaining question is therefore, Can the spectroscopy be developed to unravel specific threatening materials. THz spectroscopy addresses the energy range of 10−130 cm−1 (1 THz = 30 cm−1). This is an intermediate energy range that lies at the lower range of the vibrational frequencies of a single molecule and at the top frequency range of phonon modes of soft organic crystals. THz studies have therefore concentrated on organic crystals8 of sugars, explosives,4,9−13 biosystems, 14,15 and slow kinetic phenomena such as evaporation.16 An obstacle for remote sensing is the water absorption band.17 Recent technological progress in both CW and pulsed sources have enhanced the study of this spectral range18 and made THz radiation a suitable candidate for remote sensing.19−21 The utility of terahertz spectroscopy relies heavily on the ability to assign and interpret the spectral features.5 Typical dynamics responsible for the spectrum are long-range intramolecular motion together with soft bending and pseudorotations. The typical excitation energy in the terahertz domain ℏω is smaller than the room-temperature energy kbT by a factor of 10. These factors lead to broad lines and a congested © 2014 American Chemical Society

Received: January 1, 2014 Accepted: January 17, 2014 Published: February 10, 2014 772

dx.doi.org/10.1021/jz402801m | J. Phys. Chem. Lett. 2014, 5, 772−776

The Journal of Physical Chemistry Letters

Letter

The additional energy induced by the external field is Eμ = −μ(t)·ϵ(t). The dipole is calculated as a sum of the product between atomic coordinates with their partial charge. As a result, an additional time-dependent force acts on each atom F⃗j = ϵ⃗(t)zj, where zj is the partial charge. The atoms will respond to this force. When the frequency of the external field matches, a molecular frequency energy will be absorbed into this mode. For small molecules, where the modes are isolated, energy will accumulate quadratically as a function of time.30 This effect can be employed to determine the mode frequencies by scanning the driving force frequency. For larger systems, this procedure is inadequate because of the increasing importance of anharmonic coupling, which invalidates the normal-mode picture. Nevertheless, by scanning the driving frequency, the energy absorbed by specific modes will partition to other modes, with the final result of a total increase in the energy of the molecular total system. The energy increase will be proportional to the incident frequency and the driving amplitude. The procedure thus adopted was to apply a Gaussian pulse at time t0

based on instantaneous bond orders between neighboring atoms.29 This force field assigns partial charges to atoms, which automatically generates the dipole tensor. In linear response, the absorption can be obtained from the autocorrelation function of the dipole tensor. For small molecules using the autocorrelation function, we obtained reasonable results depending on the quality of the force field. For large systems, the quality of the spectrum degraded and became meaningless. We therefore developed a direct method for calculating the THz absorption spectrum. The method was tested on slabs of cyclotrimethylenetrinitramine (RDX) and triacetone triperoxide (TATP) crystals and was compared to experiment. The strong applied field allows one to simulate a multidimensional spectroscopy, in this case, the multispectral response of the induced dipole of the crystal in different polarization from the incident excitation. Theory: Direct Simulation of Radiation Absorption. When light propagates through matter, its oscillating electric field induce forces on the molecular dipoles. The delayed response of these dipoles back reacts to extract energy from the propagating beam. The power loss from the radiation field as a function of radiation frequency is known as the absorption spectrum. A direct simulation of this process can be carried out by adding an oscillating electric field to the MD force field.30 A pulse of radiation with the desired carrier frequency is applied to the simulated crystal. The spectrum is then obtained by evaluating the power dissipation as a function of frequency. This procedure relies on sufficient nonlinear interactions that transform the external driving field to heat, which requires a sufficient internal heat bath to absorb the input energy. For a small molecular system, this requirement will not be fulfilled. On the other hand, for a large molecular crystal, soft phonon modes create a sufficiently large internal heat bath that is able to dissipate the absorbed energy. When the radiation field is weak, linear response theory can be employed.31 Linear response relates the absorption cross section to the freely propagating dipole autocorrelation function σ(ω) =

2πω 3ℏc

2

ϵ(t ) = A e−(t − t0)

/2σ 2

cos(ωt )

where A is the pulse amplitude, σ the fwhm, and ω the carrier frequency. Repeating the calculation with varying carrier frequency and accounting for the energy added to the system as a function of frequency leads to the absorption spectrum. The amplitude A is controlled to be in the linear regime but still strong enough to overcome the random noise in the simulation. Too strong an amplitude can cause decomposition of the molecules studied. Figure 1 shows excitation by a 5 ps pulse with an accompanying energy increase in the molecular slab of RDX. The spectral resolution will depend on the pulse duration. A 5 ps pulse will lead to a resolution of ∼0.6 cm−1. Two-Dimensional (2D) THz Spectroscopy. An important consequence of the propagating electric field is the dynamical response of the charge distribution of the molecular crystal. In

∞

∫−∞ ⟨μ̂ (0)μ̂ (t )⟩eiωt dt

(1)

where μ̂ is the dipole operator and ω the absorption frequency. The input to eq 1 is the average dipole correlation function, C(t) = ⟨μ̂(0)μ̂(t)⟩. Its calculation requires a dynamical framework. Linear response together with the harmonic approximation has developed to become the primary method to simulate spectra. Linear response in a MD setting has been applied to an isolated small molecular system such as acetone, where we obtained reasonable results. The quality of the results depended primarily on the quality of the potential used. However, for a large system and for soft modes, linear response theory becomes questionable. The main reason is that classical dynamics becomes chaotic on a very short time scale, long before the dipole correlation can establish a definite frequency. For a large extended system, a direct method based on molecular dynamical simulation is almost the only practical solution. Here too, the quality of the results depends on the choice of the potential energy surface. In the present study, we used the ReaxFF force field for this task.29 The structure of this potential allows one to keep track of the partial charge on each atom. As a result, the dynamics of multiple electric moments of the matter can be evaluated at each step.

Figure 1. (Bottom) A typical excitation field ϵ(t) as a function of time with carrier frequencies of 3 (black) and 1 THz (red) The pulse duration was 5 ps. (Top) The energy accumulation as a function of time in a RDX slab of two unit cells. Black corresponds to 3 THz and red to 1 THz.

turn, the oscillating induced dipoles are a source of secondary radiation. The emitted field can be measured outside of the crystal for each input driving frequency. Most of the generated radiation is in the same frequency and polarization of the incident field. The unique part is the emitted radiation in 773

dx.doi.org/10.1021/jz402801m | J. Phys. Chem. Lett. 2014, 5, 772−776

The Journal of Physical Chemistry Letters

Letter

anharmonic character of soft organic crystals in the lowfrequency THz regime. THz Spectrum of Crystalline RDX and TATP. The theoretical methods are demonstrated by the calculation of the THz spectrum of crystalline RDX and TATP. In the case of RDX, the simulation box includes two unit cells of RDX with periodic boundary conditions. The initial coordinates were obtained from experimental data.36 The TATP crystalline structure used was isomorph 1d obtained from Reany et al.37 Two unit cells were included in the computation box. The slab was then annealed for 100 ps, reaching a temperature of 10 K. Next, the crystal was heated to 300 K using an NPT protocol. This procedure formed the initial state for which the radiation pulse was applied. The energy absorbed was calculated with an NVE protocol as the change of energy between the initial and final value. A pulse of 5 ps was employed with a varying carrier frequency ω. The integrated energy absorbed from the pulse was set to be in the linear regime on the order of 0.01 eV. Figure 2 displays a comparison between the calculated spectrim for RDX and TATP and the corresponding experimental data. The spectrum of crystalline RDX was also simulated based on a normal-mode analysis MD approach by Boyd and Boyd.39 For both solids, the low-frequency absorption spectrum can be assigned to the highest phonon modes of the molecular crystal. The peaks above 100 cm−1 correspond to the low-frequency individual torsional modes of the molecules in the slab. The frequency shifts between the calculated spectra, and the experimental data of these molecular modes are due to inaccuracies in ReaxFF. The multidimensional spectra for RDX are demonstrated in Figure 3 The incident polarization direction was chosen as the X direction. The calculated induced dipole of the molecular slab as a function of time is shown in the left panel of Figure 3. In the time domain, the first appearing dipole response is in the polarization of the incident pulse. After some delay, nonlinear response of an induced dipole appears in the perpendicular directions. In the frequency domain, the right panel in Figure 3, the diagonal signal dominates, resembling the direct absorption spectrum of Figure 2. In addition, off-diagonal peaks of up

Figure 2. The calculated THZ spectrum of crystalline RDX above and of TATP below in blue. The calculated spectrum is superimposed on the experimental spectrum for RDX4 and TATP of Wilkinson38 in dashed red.

different polarization and frequency from that of the incident pulse. This emission is due to nonlinear effects of mode-tomode coupling in the crystal. The time-dependent emitted radiation field is proportional to the dipole acceleration (d2μ)/ (dt2), which is calculated during the MD simulation. To obtain the frequency response, the dipole acceleration for each driving field ω is Fourier-transformed using the filter diagonalization procedure to obtain the response in the frequency domain f jk(ω,ω′).32,33 The spectrum of the emitted field for each input frequency and polarization j constitutes a 2D driving frequency versus response frequency and polarization k spectrograph f jk(ω,ω′). Experimentally, this response can be measured with good signal-to-noise in perpendicular polarization directions relative to the incident beam because only nonlinear response can lead to a signal in these polarization directions. The suggested one-pulse, 2D THz spectra is different from the well-developed three-pulse, third-order nonlinear 2D spectrum developed for the IR and visible spectroscopy.34,35 It is a strong field effect that is intimately related to the

Figure 3. The calculated 2D THz spectrum of the RDX crystal where the incident direction of propagation is along the X direction of the crystal with linear polarization along the Z direction. The emitted response is calculated in the three polarization directions. (Left) Time domain dipole response for different excitation frequencies. (Right) 2D spectrum using a filter diagonalization procedure to calculate the response in the frequency domain as a function of the excitation frequency. 774

dx.doi.org/10.1021/jz402801m | J. Phys. Chem. Lett. 2014, 5, 772−776

The Journal of Physical Chemistry Letters

Letter

Figure 4. The calculated 2D THz spectrum of the TATP crystal, where the incident direction of propagation is along the X direction of the crystal with linear polarization along the Z direction. The emitted response is calculated in the three polarization directions. (Left) Time domain dipole response for different excitation frequencies. (Right) 2D spectrum using a filter diagonalization procedure to calculate the response in the frequency domain as a function of the excitation frequency.

■

conversion appear where a low-frequency driving field leads to a high-frequency response, for example, the peak at the driving frequency of 10 cm−1 with a response at ∼50 cm−1 in the Z polarization. A down conversion peak can be observed in the X and Y polarizations around the driving frequency of 100 cm−1. These peaks can serve as a unique fingerprint of crystalline RDX. The calculated nonlinear 2D spectrum of TATP is shown in Figure 4. The spectrogram shows events of up conversion and down conversion. The nonlinear features are weaker in TATP compared to those in RDX. The reason is a more isotropic long-range force in the TATP crystal. The present Letter demonstrates the use of molecular dynamical simulations to calculate directly the absorption spectrum of solids by monitoring the dissipation of an electromagnetic external field. This method avoids the common assumption in spectroscopy of normal modes and weak field excitation. Direct evidence of mode-to-mode coupling is the 2D spectrogram, where the emission in a different polarization or direction is observed. The quality of the simulation depends on the quality of the potential used. For THz spectroscopy of soft organic crystals, the relevant component of the potential is the van der Waals part. The ReaxFF force field, which was used to simulate the spectrum of RDX and TATP, was calibrated for high-energy bond breaking and forming. Nevertheless, the agreement with experiment of the simulated spectrum is encouraging. An improved calibration for dispersion forces is expected to lead to a computational method with predictive power. Force fields with a better description of dielectric response have been developed.40,41 The ability to perform experimentally and simulate multidimensional spectroscopy should significantly enhance the screening ability of THz spectroscopy.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We wish to thank the support of the Department of Homeland Security, ALERT Center of Excellence Award 2013-ST-061ED0001.

■

REFERENCES

(1) Willer, U.; Schade, W. Anal. Bioanal. Chem. 2009, 395, 275−282. (2) Pacheco-Londoño, L. C.; Ortiz-Rivera, W.; Primera-Pedrozo, O. M.; Hernández- Rivera, S. P. Anal. Bioanal. Chem. 2009, 395, 323− 335. (3) Davies, A. G.; Burnett, A. D.; Fan, W.; Linfield, E. H.; Cunningham, J. E. Mater. Today 2008, 11, 18−26. (4) Leahy-Hoppa, M. R.; Fitch, M. J.; Osiander, R. Anal. Bioanal. Chem. 2009, 395, 247−257. (5) Kemp, M. C. IEEE Trans. Terahertz Sci. Technol. 2011, 1, 282− 292. (6) Jepsen, P. U.; Cooke, D. G.; Koch, M. Laser Photonics Rev. 2011, 5, 124−166. (7) Hoffmann, M. C.; Yeh, K.-L.; Hebling, J.; Nelson, K. A. Opt. Express 2007, 15, 11706−11713. (8) Motley, T. L.; Korter, T. M. Chem. Phys. Lett. 2008, 464, 171− 176. (9) Allis, D. G.; Korter, T. M. ChemPhysChem 2006, 7, 2398−2408. (10) Whitley, V. H.; Hooks, D. E.; Ramos, K. J.; OHara, J. F.; Azad, A.; Taylor, A.; Barber, J.; Averitt, R. Anal. Bioanal. Chem. 2009, 395, 315−322. (11) Whitley, V. H.; Hooks, D. E.; Ramos, K. J.; Pierce, T. H.; OHara, J. F.; Azad, A. K.; Taylor, A. J.; Barber, J.; Averitt, R. D. J. Phys. Chem. A 2011, 115, 439−442. (12) Barber, J.; Hooks, D. E.; Funk, D. J.; Averitt, R. D.; Taylor, A. J.; Babikov, D. J. Phys. Chem. A 2005, 109, 3501−3505. (13) Wilkinson, J.; Konek, C. T.; Moran, J. S.; Witko, E. M.; Korter, T. M. Chem. Phys. Lett. 2009, 478, 172. (14) Korter, T.; Plusquellic, D. F. Chem. Phys. Lett. 2004, 385, 45− 51. (15) Plusquellic, D. F.; Siegrist, K.; Heilweil, E. J.; Esenturk, O. ChemPhysChem 2007, 8, 2412−2431.

ASSOCIATED CONTENT

S Supporting Information *

RDX and TATP unit cell information and ReaxFF-lg force field parameters. This material is available free of charge via the Internet at http://pubs.acs.org. 775

dx.doi.org/10.1021/jz402801m | J. Phys. Chem. Lett. 2014, 5, 772−776

The Journal of Physical Chemistry Letters

Letter

(16) Innocenzi, P.; Malfatti, L.; Piccinini, M.; Sali, D.; Schade, U.; Marcelli, A. J. Phys. Chem. A 2009, 113, 9418−9423. (17) Withayachumnankul, W.; Fischer, B. M.; Abbott, D. Proc. R. Soc. London, Ser. A 2008, 464, 2435−2456. (18) Ferguson, B.; Zhang, X.-C. Nat. Mater. 2002, 1, 26−33. (19) Mittleman, D.; Gupta, M.; Neelamani, R.; Baraniuk, R.; Rudd, J.; Koch, M. Appl. Phys. B: Laser Opt. 1999, 68, 1085−1094. (20) Williams, B. S.; Kumar, S.; Hu, Q.; Reno, J. L. Opt. Express 2005, 13, 3331−3339. (21) Padilla, W. J.; Taylor, A. J.; Highstrete, C.; Lee, M.; Averitt, R. D. Phys. Rev. Lett. 2006, 96, 107401. (22) Oxley, J.; Smith, J.; Brady, J.; Dubnikova, F.; Kosloff, R.; Zeiri, L.; Zeiri, Y. Appl. Spectrosc. 2008, 62, 906−915. (23) Brauer, B.; Dubnikova, F.; Zeiri, Y.; Kosloff, R.; Gerber, R. B. Spectrochim. Acta, Part A 2008, 71, 1438−1445. (24) Allis, D. G.; Prokhorova, D. A.; Korter, T. M. J. Phys. Chem. A 2006, 110, 1951−1959. (25) Korter, T.; Balu, R.; Campbell, M.; Beard, M.; Gregurick, S.; Heilweil, E. Chem. Phys. Lett. 2006, 418, 65−70. (26) Hakey, P. M.; Allis, D. G.; Ouellette, W.; Korter, T. M. J. Phys. Chem. A 2009, 113, 5119. (27) Day, G. M.; Zeitler, J.; Jones, W.; Rades, T.; Taday, P. J. Phys. Chem. B 2006, 110, 447−456. (28) Saito, S.; Ohmine, I. J. Chem. Phys. 1998, 108, 240−251. (29) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396−9409. (30) Bowman, J. M.; Zhang, X.; Brown, A. J. Chem. Phys. 2003, 119, 646. (31) Kubo, R. J. Phys. Soc. Jpn. 1957, 12, 570−586. (32) Pang, J. W.; Dieckmann, T.; Feigon, J.; Neuhauser, D. J. Chem. Phys. 1998, 108, 8360. (33) Mandelshtam, V. A. Prog. Nucl. Magn. Reson. Spectrosc. 2001, 38, 159−196. (34) Tanimura, Y.; Mukamel, S. J. Chem. Phys. 1993, 99, 9496. (35) Zanni, M. T.; Asplund, M. C.; Hochstrasser, R. M. J. Chem. Phys. 2001, 114, 4579. (36) Choi, C.; Prince, E. Acta Crystallogr., Sect. B 1972, 26, 2857− 2862. (37) O. Reany, M. B.; Kzpon, M.; Keinan, E. Cryst. Growth Des. 2009, 9, 3661−3670. (38) Wilkinson, J.; Konek, C. T.; Moran, J. S.; Witko, E. M.; Korter, T. M. Chem. Phys. Lett. 2009, 478, 172−174. (39) Boyd, S. G.; Boyd, K. J. J. Chem. Phys. 2008, 129, 134502. (40) Karasawa, N.; Goddard, W. A., III. Macromolecules 1995, 28, 6765−6772. (41) Tkatchenko, A.; Scheffler, M. Phys. Rev. Lett. 2009, 102, 073005.

776

dx.doi.org/10.1021/jz402801m | J. Phys. Chem. Lett. 2014, 5, 772−776