**UFRJ**

# Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain

(Parte **1** de 4)

Molecular Dynamics Simulation of Heat Conduction through a Molecular Chain

Christian Schroder,*,† Vyacheslav Vikhrenko,‡ and Dirk Schwarzer§

Instituteof ComputationalBiologicalChemistry,UniVersity of Vienna, Vienna, Austria,BelarusianState TechnologicalUniVersity, Minsk, Belarus,and Max-Planck-Institutefor BiophysicalChemistry,Gottingen,Germany

ReceiVed: April 17, 2009; ReVised Manuscript ReceiVed: September 14, 2009

This work deals with a moleculardynamicssimulationanalysisof the intramolecularvibrationalenergy transfer in a system of two chromophores, azulene and anthracene, bridged by an aliphatic chain and is motivated by corresponding laser experiments. After selective excitation of the azulene chromophore, the subsequent intramolecular vibrational energy redistribution is monitored by analyzing the transient temperatures of the two chromophores and the chain between them. The main focus concerns the heat conduction process in the chain. Therefore, the chain length was varied from 0 to 19 CH2 units. In addition, methoxymethyl, 1,2- dimethoxyethyl, and a thiomethoxymethyl chains were studied. The investigation of the intramolecular vibrational energy process was decomposed into a temporal analysis and a spatial analysis. For short alkyl chains, the time constant of energy relaxation increases proportionally to the chain length. However, for longer chains, the time constant characterizing the energy decay of the azulene chromophore saturates and becomes independent of the chain length. This behavior is consistent with experimental findings. The spatial analysisshows more or less exponentialdecay of the temperaturealong the chain near the excitedchromophore. In additional simulations, the two chromophores were thermostatted at different temperatures to establish a constant heat flux from the azulene to the anthracene side. The steady-state temperature profiles for longer alkyl chains show strong gradients near the two chromophores and constant but weak gradients in the central part of the chain. Both simulation methods indicate that strong Kapitza effects at the boundaries between each chromophore and the molecular chain dominate the intramolecular energy flux.

Introduction

Both intramolecular vibrational energy redistribution (IVR) and intermolecular vibrational energy transfer (VET) are of fundamental importance in gas- and condensed-phase chemical dynamics. These nonreactive elementary processes determine the rates, pathways, and efficiencies of chemical reactions. Whenever a polyatomic molecule is locally excited by some chemical or electromagnetic activation process, IVR provides the basic mechanism of energy exchange between different vibrational degrees of freedom (i.e., normal modes). Here, the term “local excitation” applies to the preparation of some zeroorder eigenstates or superpositions thereof or even to a local distortion of the molecule. In the presence of a solvent bath, excess energy is transferred irreversibly to the bath degrees of freedom.1,2 In addition, the solvent might affect the intramolecular redistribution of energy, which is usually taken into account by the notion of “solvent-assisted” IVR.1,3-6 Whereas the assumption of rapid IVR (on the time scale of reactive processes) is at the heart of statistical unimolecular rate theories,7,8 intermolecular equilibration determines the applicability of canonical transition state theory.9 Often, a separation of time scales for IVR and VET is assumed, such as in Rice-Ramsberger-Kassel-Marcus Theory.7 In the gas phase, where the average time between collisions is rather long, this assumption is usually correct, at least for low densities.10,1 In solution, however, it is questionable in general.1,5,12 Because IVR is driven by state-specific anharmonic couplings and because dissipative intermolecular energy transfer is a strongly frequency-dependent phenomenon, favoring VET by lowfrequencyvibrations,IVR might be incompleteon the time scale of energy-transferring collisions. In fact, in recent molecular dynamics simulations, this phenomenon was demonstrated for

CH2I2 in CDCl3 solvent.12 For large molecules, where, because of the high density of states, IVR is often no longer strictly state-specific but statistical, the separability of time scales can usually safely be assumed even in liquid solution.3,13,14

Traditionally, IVR is discussed by reference to the eigenstate space of some zero-ordermolecularHamiltonian.15-17 The study of IVR has been extended to the investigation of vibrational energy transport between different chromophore parts of a large molecule through a flexible molecular chain.4,18,19,20 In such a case, in addition to the energy and normal-mode picture, a (Cartesian)spatial dimensionof IVR comes into play. In a series of experiments,4,18 two chromophores are attached to an alkyl chain as illustrated in Figure 1. The first optical chromophore, azulene, is excited by a femtosecond laser pulse into its first excited singlet state, S1, from which it returns via internal conversion within e1p st ot he S0 state,21,2 creating a vibrationally highly excited ground-state azulene moiety. Subse- quently, the excess energy flows through the alkyl chain to the second chromophore (e.g., benzene, 1-methylene anthracene, or anthracene).23 The overall process can be monitored via transient absorption of photons from a second laser pulse tuned to the red edge of the azulene or the second chromophore (e.g., anthracene) spectra. From the change in the optical density, ∆OD, of the chromophores, the temperatures of the azulene and anthracene can be calculated.4,23 In addition to IVR, the intramolecularenergy flow is also characterized by a “harmonic dephasing” process of the eigenfrequencies of the used com-

* To whom correspondence should be addressed. E-mail: christian@ mdy.univie.ac.at.

† University of Vienna. ‡ Belarusian State Technological University. § Max-Planck-Institute for Biophysical Chemistry.

10.1021/jp903546h 2009 American Chemical Society Published on Web 1/24/2009

pound. In a strict sense, this process is not IVR because, theoretically,no energy is transferredbetweenthe normal modes of the compound. However, because of the time scale of the experiment and the anharmonicities, the initial excited state is not restored. As a result, the harmonic dephasing (which occurs on a subpicosecondtime scale) contributesto the energy transfer between the two chromophores within the experimentally accessible observation time. The IVR time constant is in the range of 1-5 ps depending on the length of the chain and the second chromophore,23 whereas VET occurs with a time constant of approximately 20-50 ps in the molecular solvent used.4,23 Because of the observed separation of time scales between IVR and VET, it can be assumed that IVR is complete on the time scale of VET. However,the VET time is remarkably close to the energy relaxation time of bare azulene in the same solvent,21 which suggests that low-frequency vibrations of the bridged compounds are not significantly populated.23

Initial experimental data on azulene-anthracene compounds with short molecular bridges showed a linear dependence of the IVR time constant on the alkyl chain length.4 These findings were confirmed by molecular dynamics simulations in which the interchromophore IVR time was found to be proportional to the alkyl chain length for up to four CH2 units.18 In principle, the linear dependence of the IVR time constant on the chain length can be explained by the Fourier law of heat conduction.24 The two chromophores can be regarded as heat baths, which, however, because of their finite heat capacities, change their temperatures as a result of energy exchange with the bridge degreesof freedom.For a homogeneousone-dimensionalspatial grid with end points that are held at different temperatures, Fourier’s heat conduction equation yields a stationary solution with a linear temperature profile and, thus, a constant flux at all grid points. If the heat reservoir (end-point) temperatures vary in time but the reservoir heat capacities are still large as compared to that of the one-dimensional bridge, a similar nonequilibrium stationary state with an approximately linear temperature field is established, with a slope that decreases in the course of temperature equilibration.4 As a result, Fourier’s law predicts temperature equilibration times proportional to the length of the intervening spatial grid. Recent results, however, suggest a saturation of the interchromophore energy-transfer time for greater bridge lengths,18 pointing to a possible breakdown of Fourier’s law of heat conduction.

An “anomalous” thermal transport has also emerged from a large number of simulation studies on heat conduction in onedimensional chains.25,26 A major problem with these studies, however, is that the results strongly depend on the methods and the boundary conditions used.

Another type of anomalous thermal transport arises from the so-called Kapitza effect and has been examined experimentally and theoretically.27-34 This effect appears at interfaces when system parameters are changing considerably, for example, at solid-liquid interfaces or at the contact of two crystallites in polycrystals.

In macroscopic theory, heat transfer across interfaces is taken into account by choosing proper boundary conditions,35,36 such as by including a thermal resistance of the interface. In a microscopicdescription,the Kapitzaeffectis consideredthrough a coordinate-dependent thermal conductivity that is larger in the bulk than in the boundary layers.

The aim of the present work is to investigate the heat conduction through an aliphatic chain between an azulene and an anthracene chromophore by means of nonequilibrium molecular dynamics simulations. This process is determined by anharmonic couplings between the vibrational degrees of freedom in the molecules. The electronic contribution to the heat transport can safely be ignored.

Two approaches are usually used for MD calculations of transport coefficients. One of them is based on linear response theoryand the Green-Kubo relationsthat representthe transport coefficients through time integrals of the corresponding time correlation functions.31,37,38 This approach is appropriate for spatially homogeneous systems because averaging over the system volume increases statistics by orders of magnitude. For spatially inhomogeneous systems, it is usually prohibitive because of problems with statistics.

In the other approach, a nonequilibrium system state is created, and transport coefficients are calculated as constants of proportionality between fluxes and gradients of the corresponding (thermodynamic) variables. In this approach, either time evolution of the system is considered, or a nonequilibrium steadystate is created.Just nonequilibriumMD simulationswere chosen in this work to investigate the energy-transfer problem.

This articleis organizedas follows:The next sectiondescribes the details of the simulations used to mimic the experimental situation. To reproduce the experimental data18,23 for the azulene-anthracene system as closely as possible, a new molecular dynamics force field (presented in Appendix A) is introduced for these bridged compounds, which fits the experimentalspectralinformationvery well.The emphasisof this work lies on the spatial characterization of the IVR, that is, the IVR time constants of the two chromophores and the temperature profile within the alkyl chain. After a short introduction to the general concept of temperature and Fourier heat conduction, the necessity of introducing spatially dependent heat conductivity coefficients due to the “Kapitza” effect is explained. The presentation of the results is subdivided into two parts: First, the transient temperature relaxation is discussed, including the dependence on the chain length and chain composition. Ad- ditionally,one simulationseries of azulene-(CH2)3-anthracene in Xe is presented. The spatial analysis deals with the spatial temperatureprofile in the case of the nonequilibriumsimulations (emulating the laser experiments) and steady-state simulations with thermostatted chromophores. Finally, we summarize the results.

Simulation Details

In the present work, intramolecular vibrational energy redistribution in a model series of azulene-(CH2)n-anthracene compounds (cf. Figure 1) with varying alkyl chain lengths (0 e n e 19) is simulated via nonequilibrium classical molecular dynamics, mimicking the experimental pump/probe investiga-

Figure 1. Definition of atom types in the azulene (azu) and anthracene (ant) moieties, as well as in the alkyl chain (alk). Hydrogen atoms are assignedthe same Greek indicesas the carbonsto which they are bound.

14040 J. Phys. Chem. A, Vol. 113, No. 51, 2009 Schroder et al.

tions performed on the same type of model systems.4,18 In accordancewith the interpretationof the experiments,the model systems are assumed to be equilibrated at a given heat bath temperature of the sample (usually 300 K, unless explicitly stated otherwise) prior to excitation by a subpicosecond visible light pulse. The latter is used in the experimentsto locally excite the azulene chromophore into its lowest excited singlet state,

S1, from which the azulene unit returns to its S0 state within e1 ps by internal conversion. Because local energy redistribu- tion in the azulene chromophore cannot experimentally be separated from the internal conversion process, it is assumed thatson the time scale of through-bridge azulene-to-anthracene energy transfersthe initial vibrational excess energy is randomized among the azulene vibrations. Because of the relatively large number of local chromophore vibrational degrees of freedom, this local energy equilibrium (which is, in principle, microcanonical, with additional dispersive contributions due to canonical equilibrium prior to excitation) can be interpreted in terms of an initial nonequilibrium temperature applying to the local azulene degrees of freedom. Therefore, the initial mean excessenergyof 2.23 eV (18000cm-1) is statisticallydistributed among the azulene vibrations according to a local canonical equilibrium. When vibrational degrees of freedom are treated classically, this heats the azulene to a temperature of approximately 850 K. Shortly thereafter, a coherent dephasing process transfers energy to the rest of the molecule on a subpicosecond time scale.39 The delocalized energy heats the alkyl chain and anthracene by approximately 20-80 K and decreases the azulene temperature to 750 K.

Initialization of configurations in phase space was thus performed as follows: Starting from the minimum-energy geometries of model systems, the azulene was thermostatted at 750 K and the rest of the molecule at 320 K for a simulation period of 100 ps using a Nose-Hoover algorithm of CHARMM.40 This thermostatting avoids the “flying icecube” problem.41 In other words, the thermostat does not introduce a translationaldrift. After this thermostattedphase, all thermostats were switched off. This ensured that there were no artificial Kapitza effects caused by the link between thermostatted and nonthermostated atoms in the system. The trajectory was propagated for an additional 150 ps using a time step of 0.1 fs. For each alkyl bridge length, an ensemble of 2000 such nonequilibriumtrajectorieswas producedto ensure the statistical quality of the data.

The time scales of intramolecular energy redistribution will certainly depend on the shape of the potential energy surface. For large molecular systems as investigated in the present work, with more than 120-300 degrees of freedom, ab initio information on the energy surface even at low energies is beyond reach. Although simulation with ab initio forces calculated “on the fly” using efficient electronic structure methods have become increasingly popular,42,43 the typical time scales covered are still in the lower picosecond range, which is too short for slow IVR. Therefore, we chose to rely on a local analytical force field, as traditionally employed in classical molecular dynamics simulations, parametrized to reproduce, as closely as possible, the spectrumof vibrationalfrequenciesobtainedfrom a combination of densityfunctionaltheory(DFT) calculationsand experimental information. As a result, the anharmonicity of the potential energy function and thus the detailed mode-to-mode couplings are, to some extent, “accidental”. However, the usual functional forms of bond, angle, and torsional potentials are realistic enough to obtain a physically reasonable parametrization of the potential energy surface. Experience from previous simulation projects on IVR and VET in solution,12 in particular the time scales obtained from molecular dynamics simulations as compared to experimental values, seems to justify our approach. The parametrization of the chromophores azulene (denoted as azu in the equations and tables) and anthracene (ant) and the alkyl chains (alk) is described in detail in Appendix A.

Temperature and Heat Conduction

In this work, many nonequilibrium MD simulations of azulene-anthracene systems were performed in order to juxtaposecomputationalresultsand experimentaldata.4,18,39Because we want to concentrateon IVR withinthe moleculein this work, most of the simulation results deal with the pure azu-alk-ant derivatives without a solvent. Consequently, VET processes are neglectedand not discussedhere. The originalexperimentaldata consist of absorption-time profiles of the azulene and anthracene chromophores as a function of the length of the alkyl chain that bridges the two chromophores,as illustratedin Figure 1. This measured change in optical density, ∆OD(t), was been related to the temperature relaxation of the chromophores.23 An assignment of individual temperatures Ti(t) to each subunit i ∈ {azu, ant} assumes that IVR in the chromophoresis much faster than the energyredistributionbetweenthem throughthe aliphatic chain.

Ti(t) is calculated in our simulations by summing the kinetic energies of the corresponding atoms j belonging to that subunit25,4 where ni is the number of atoms in the corresponding subunit. In accordance with the experiments, the chromophore temper- atures Tazu(t) and Tant(t) in the simulations relax, to a good approximation monoexponentially, to an equilibrium tempera- ture Teq

The time constantsτi characterizethe energyequilibrationwithin the molecule and can be directly compared with experiment.

After initial excitation, the local temperature of the azulene is is negative,because anthracenegains energy until Teq is reached. To elucidate the mechanism of heat transfer in detail, spatial temperature profiles through the molecular chain were calcu- for each CH2 unit k,0 e k e n, was analyzed separately (l is the carbon-carbon bond length in the chain).

Based on a classic energy conservation equation for a chain in continuum representation

the heat flux J(x,t) distributes the local excess temperature until an equilibrium is achieved. In this equation, cL is the specific heat per unit length. The classical way to describe the heat conduction is according to Fourier’s law of heat conduction,24

3nikB ∑ j

Heat Conduction through a Molecular Chain by MD J. Phys. Chem. A, Vol. 113, No. 51, 2009 14041 which states that the heat flux, J(x,t), is proportional to the temperature gradient

The heat conduction coefficient κ and its inverse, the heat resistance r ) 1/κ, are usually considered as constant. In this case, eq 5 results in a diffusive process of heat transfer with a

uniform thermal diffusivity R) κ/cL. However, the molecular systemstreatedhere are characterizedby spatialinhomogeneities and hence require spatially varying heat resistances r ) r(x) system, the heat capacity can be approximated by cL ) 9kBl-1 .

The heat flux can be determined by the power acting from one part of the molecule to another, i.e.:

Here, vj is the velocity of the jth atom, and Fj is the sum of all forces acting on atom j due to the intramolecular potentials

(bond, angle, dihedral, improper torsions) at the corresponding interface. Equation 6 can be used to calculate the heat flux at any bond in the alkyl chain.

Results and Discussion

We present our results in three sections. The first covers the transient temperature relaxation of the two chromophores after the initial excitation of azulene. This part closely resembles the experimental situation, and the corresponding time constants can be directly compared. The second part is concerned with the temperature profiles of the chain during relaxation. Finally, we discuss the results of simulations with chromophores thermostatted at different temperatures to establish a steadystate energy flux.

Time-Variant Temperature Relaxation. Because the attached optical chromophores azulene and anthracene are relatively large and the energy flux through the chain is relatively slow, local microcanonical equilibrium among their respective vibrational degrees of freedom can be described by an effective temperature,as in eq 1. The temperaturewas calculatedfor each subunit at each simulation step and averaged over the corresponding ensemble of trajectories. We start our analysis with nonequilibrium simulations of azu-(CH2)3-ant in a bath of 225 Xe atoms. After initial excitation of the azulene, the energy was redistributed among all atoms of the azu-(CH2)n-ant compoundvia IVR and finally transferred to the heat bath. The corresponding temperature relaxation Ti(t) is depicted in Figure 2. Likewise, we also calculatedthe temperatureof the alkylchain.However,as shown below, during relaxation, the vibrational degrees of freedom of the chain were not at all equilibrated, which makes the assignment of a temperature questionable. Rather, it can be understood as a measure of the internal energy of the chain.

The temperaturesof azulene and anthracenecan be well fitted by means of biexponential functions with equilibrium temperatures of 300 K. The longer time constant of approximately τVET ) 100 ps is attributed to the transfer of vibrational energy to the bath (VET) and is in reasonable agreement with the experimental τVET value of 35 ps.23 The faster time constants, τIVR, attributed to the intramolecular energy flow turn out to be very similar to the values calculated for the isolated molecule without solvent; that is, the environment does not seem to influence the intramolecular process significantly. This is consistent with the experimental observation where no depen- dence of τIVR on the solvent was found.23 Therefore, most of the calculations were done for the isolated molecules.

The results for the isolated azu-(CH2)n-ant model systems with n ) 0-19 methylene units are presented in Figure 3.

During IVR, the excess energy is transferred from the azulene to the rest of the molecule. Consequently, the temperature of the azulene decreases, and the temperatures of the chain and anthracene increase until Teq is reached. Because the excess energy is redistributed solely among the atoms of the azu-(CH2)n-ant unit and not transferred to a heat bath, the equilibrium temperature Teq is significantly higher than 300 K (cf. Figure 3). It decreases with chain length because of the increase in the heat capacity of the alkyl chain. For n g 2, the temperature dynamics of each azu-(CH2)n-ant compound in Figure 3b-f can be fitted monoexponentially according to eq

3. Molecules with directly linked chromophores or with just one methylene group between them (n ) 0, 1) require a biexponential fit to model the energy equilibration. This is clearly visible in the case of azu-ant in Figure 3a. The time constants are summarized in Table 1. From the viewpoint of the intramolecularMD potential of these molecules, the azulene atoms are directly coupled to anthracene atoms by angle and dihedral potentials for n < 2. The number of potential terms affecting both chromophores is larger for n ) 0 than for n ) 1.

(Parte **1** de 4)