J. Phys. Chem. B 2008, 112, 313-320
313
Umbrella Sampling of Solute Vibrational Line Shifts in Mixed Quantum-Classical Molecular Dynamics Simulations† Christine M. Morales and Ward H. Thompson* Department of Chemistry, UniVersity of Kansas, Lawrence, Kansas 66045 ReceiVed: June 28, 2007; In Final Form: October 15, 2007
An umbrella sampling approach for vibrational frequency line shifts is presented. The technique allows for efficient sampling of the solvent configurations corresponding to frequency shifts of a solute in mixed quantumclassical simulations. The approach is generally applicable and can also be used within traditional perturbation theory calculations of frequency shifts. It is particularly useful in the extraction of detailed mechanistic information about the solute-solvent interactions giving rise to the frequency shifts. The method is illustrated by application to the simple I2 in a liquid Xe system, and the advantages are discussed.
1. Introduction Molecular dynamics (MD) simulations are a powerful tool for modeling vibrational spectroscopies of molecules in condensed phase environments. Not only can the vibrational frequency shifts and line broadening due to fluctuations in the surroundings be accurately described, a detailed molecular-level understanding of the underlying mechanisms can also be obtained. This is because MD allows one to correlate, for example, the instantaneous frequency shift with the configuration of the surrounding solvent molecules. Moreover, the contribution to the frequency shift of each solvent molecule can be calculated.1-3 We have recently demonstrated the level of mechanistic information that can be obtained.3 Specifically, the mechanism(s) of solvent-induced vibrational frequency shifts can be determined as a function of the shift. Thus, common models for describing frequency shifts, such as the isolated binary collision model,5-7 can be directly and stringently tested by determining how many solvent molecules contribute to the instantaneous line shift. The groundwork can also be laid for the development of better models by evaluating where the largest-contributing solvent molecules are with respect to the solute molecule. While the simulation of a linear infrared or Raman spectra for comparison with experimental measurements does not generally require enhanced sampling, obtaining this kind of detailed mechanistic information can require greater challenges. Specifically, mechanistic details often stand out most prominently when the instantaneous vibrational line shifts are large (in absolute value), and it can be useful to understand whether models, for example, the isolated binary collision model, are accurate in the limit of large shifts. Thus, a full understanding of the origins of the line shifts can require sampling of rarely achieved frequency shifts (vide infra). It is worth noting, however, that this is not the only motivation for methods providing enhanced sampling of vibrational line shifts. Sufficient sampling of solute frequency shifts may be difficult to achieve in heterogeneous systems such as interfaces or nanoconfined †
Part of the “James T. (Casey) Hynes Festschrift”. * To whom correspondence should be addressed. E-mail: wthompson@ ku.edu.
solvents, where the solute frequency shifts can arise from several different, long-lived local environments. An analogous demand also arises when modeling a solute that can undergo a chemical reaction or conformational change. In those cases, umbrella sampling of the frequency shift may facilitate determination of the full spectrum resulting from reactants, products, and transition states, and the line shift may, in some cases, serve as a proxy for a reaction coordinate. Umbrella sampling8,9 is a means to generate trajectories with sampling over an extended range of line shifts or targeted sampling over a specific narrow range. Here, we show that umbrella sampling efficiently accesses mechanistically relevant regions of configuration space and is relatively simple to implement within the context of a mixed quantum-classical simulation of vibrational line shifts. While we apply the umbrella sampling within the context of a mixed quantum-classical MD approach, wherein the solute vibrational wave function influences the solvent dynamics and vice versa, as discussed in the Appendix, it may equally well be used with the traditional perturbation theory approach,10-12 for which the solute vibrational mode is fixed during the MD simulation. In section 2, we describe the implementation of umbrella sampling based on the vibrational frequency in the context of mixed quantum-classical molecular dynamics; the formulation of the method for standard perturbation-theory-based molecular dynamics simulations is given in the Appendix. An illustrative application of this technique to the analysis of I2 vibrational frequencies in liquid Xe is presented in section 3. Special attention is paid to the applicability of the umbrella sampling technique as a way to investigate the molecular-level mechanisms of vibrational frequency shifts in condensed phases. Concluding remarks are offered in section 4. 2. Methods 2.1. Mixed Quantum-Classical Molecular Dynamics. We have previously used vibrationally adiabatic, mixed quantumclassical molecular dynamics to analyze the microscopic origins of frequency shifts2,3 and vibrational energy transfer.9,13 This approach, which serves as the basis of the umbrella sampling technique presented in this paper, has been described in detail elsewhere.14 The mixed quantum-classical approach includes feedback between the solute vibration and solvent dynamics,
10.1021/jp075038d CCC: $40.75 © 2008 American Chemical Society Published on Web 12/15/2007
314 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Morales and Thompson
which increases when the instantaneous frequency shift is large. Such an approach may be particularly useful when umbrella sampling is used to bias the trajectory toward large frequency shifts. We also note that an analogous frequency umbrella sampling treatment is possible under a standard perturbation theory approach, as outlined in the Appendix. Briefly, the full classical Hamiltonian for the I2 in Xe system can be written as
H(r, pr, e, pe, Q, P) )
p2r
+
2µ
p2e
N
+
2µr2e
P2j
+ V(r, e, Q) ∑ j)2 2m
(1)
j
where r and e are the bond distance and orientational coordinates of the diatomic solute with corresponding momenta pr and pe, Q ) (Q1, Q2, ..., QN) are the coordinates of the solute centerof-mass and N - 1 solvent atoms with corresponding momenta P ) (P1, P2, ..., PN), µ is the solute reduced mass, and mj is the mass of solvent atom j. The potential energy, V(r, e, Q), describes the solute vibrational potential, solute-solvent interactions, and solvent-solvent interactions; details are given in section 3.1. Using a mixed quantum-classical adiabatic approximation for the I2 vibration, the “fast” solute vibration is described quantum-mechanically, while all other degrees of freedom are governed by classical equations of motion. The quantum mechanical Hamiltonian hˆ r for a given set of center-of-mass positions, e and Q, follows accordingly
hˆ r(e, Q) )
pˆ 2r 2µ
+ Vˆ (rˆ; e, Q)
(2)
The vibrational Schro¨dinger equation is solved at each time step for the adiabatic vibrational states, φn and corresponding energies, En
hˆ r(e, Q)φn(r; e, Q) ) En(e, Q)φn(r; e, Q)
(3)
These solutions depend implicitly on time through the changes in the classical variables e and Q. A potential-optimized15,16 discrete variable representation17-20 (DVR) basis set consisting of 30 grid points is used for all quantum mechanical calculations. Only the three lowest vibrational energy levels are considered, and the Lanczos21,22 method is used to diagonalize the vibrational Hamiltonian at each time step. A vibrationally adiabatic Hamiltonian is then generated with the solved eigenfunction and energy level for state n N
Hn(e, pe, Q, P) )
P2j
+ ∑ j)2 2m
2µr2e
P2j
p2e
j
N
)
p2e
+ ∑ j)2 2m j
2µr2e
+ 〈φn|hˆ r|φn〉r
sented here can also be used without this approximation with some additional algorithmic and computational effort. The largest effect of this approximation is a reduction in the calculated magnitude of vibrational red shifts;3 mechanistic analysis hinges on contributions to hˆ r from solute-solvent interactions, which should not be substantially altered. We have shown previously that the vibrational line shifts at each time step can be decomposed into contributions from particular solute-solvent interactions to provide detailed insight into the mechanism of the shifts.3 Relative to the gas-phase frequency, the instantaneous fundamental vibrational line shifts ∆ω are directly obtained from the energy gap ∆E01 at each time step. They are decomposed into pairwise additive solute-solvent potential terms ∆Vj01/p and nonadditive terms that depend implicitly on e and Q and are grouped together as ∆ω(0) gas ∆ω ) [∆Esolv 01 - ∆E01 ]/p
[∑ ] N
∆ω ) ∆ω(0) +
∆V(j) 01 /p
j)2
where
∆ω(0) ) [{∆Tr01 - ∆Tr01(gas)} + {∆Vˆ AB 01 (r) ∆Vˆ AB 01 (r)(gas)}]/p (6) Analysis of the additive terms ∆Vj01/p results in detailed information about the solvent configurations responsible for the distribution of instantaneous line shifts. However, significant sampling is required, which is not easily accomplished for the more rarely achieved shiftsstypically those that are large in magnitudesor if information is desired about a narrow frequency range. 2.2. Umbrella Sampling for Vibrational Line Shifts. Here, we describe the implementation of umbrella sampling8,23 to sample the configuration space around a given vibrational line shift. We introduce a biasing potential ∆V(i) umb into the Hamiltonian to selectively sample configurations with line shifts near a chosen value, ∆ω(i) umb. This biasing term is harmonic in form
1 2 (i) ∆V(i) umb ) k[p(∆ω - ∆ωumb)] 2
(7)
A Lagrange multiplier, λ, has also been introduced to ensure that the directional vector e remains normalized. The equations of motion are thus determined by the following modified Hamiltonian (i) H ˜ (i) n (e, pe, Q, P) ) Hn(e, pe, Q, P) + ∆Vumb + λ(e‚e - 1) N
+ En(e, Q)
(5)
(4)
which determines the equations of motion for the classical variables. In the present work, the classical trajectory is propagated in the ground state, n ) 0. We emphasize that, in contrast to perturbation theory methods, this approach allows the solute vibrational wave function to influence the solvent dynamics just as the solvent dynamics influence the solute vibrational wave function. We also note that in order to simplify the molecular dynamics integration algorithm, we have invoked a rigid rotor approximation for the solute, in which p2e /(2µrˆ2) is replaced by p2e /(2µr2e ). The umbrella sampling approach pre-
)
P2j
+ ∑ j)2 2m j
1
p2e 2µr2e
+ En(e, Q) +
2 k[p(∆ω - ∆ω(i) umb)] + λ(e‚e - 1) (8)
2
On the basis of the modified Hamiltonian, forces are evaluated according to Hamilton’s equations of motion. For a groundstate (n ) 0) trajectory sampling the fundamental vibrational frequency, this gives
e3 )
pe ∂H ˜ (i) 0 ) 2 - 2λe ∂pe µr
(9)
Umbrella Sampling of Solute Vibrational Line Shifts
Q4 j )
J. Phys. Chem. B, Vol. 112, No. 2, 2008 315 TABLE 1: Potential Parameters
∂H ˜ (i) Pj 0 ) ∂Pj mj
(10)
[
]
[
]
Lennard-Jones23
∂H ˜ (i) ∂E1 ∂E0 ∂E0 0 (∆E01 - p∆ω(i) p3 e ) )-k umb) ∂e ∂e ∂e ∂e (11) ∂H ˜ (i) ∂E0 ∂E1 ∂E0 0 )-k (∆E01 - p∆ω(i) P4 j ) umb) ∂Qj ∂Qj ∂Qj ∂Qj (12) We note that the additional term generated by the umbrella potential in the equation for p3 e (P4 j) depends on the gradient (∂En/∂e) ((∂En/∂Qj)), which is obtained using the HellmannFeynman theorem, for example, (∂En/∂e) ) 〈φn|(∂Vˆ /∂e)|φn〉r. In this expression, the force operator, for example, (∂Vˆ /∂e) expresses the classical force as a function of the diatomic bond length, in the DVR basis. Thus, the umbrella potential only requires the additional calculation of the gradients of E1, which requires essentially no additional computational effort. We use a simple Verlet leapfrog algorithm to integrate the equations of motion. If the centrifugal kinetic energy of the diatomic solute is to be included in hˆ r, the equations of motion must be solved self-consistently, as described previously.14 From each simulation, a frequency distribution for the biased potential can be obtained in the form of a histogram. For a given (i) ∆V(i) umb, Punbiased(∆Ω) can be considered valid within a limited “window” of ∆Ω values near ∆ω(i) umb, for which the probability distribution is accurate. To obtain the unbiased frequency distribution over a wider range, a series of simulations are carried out with different values for ∆ω(i) umb, such that the windows overlap to cover the desired range of frequencies.24 An unbiased histogram is obtained by reweighting the probability according to the biasing potential
〈e[β∆Vumb-f(i)] × δ(∆ω - ∆Ω)〉biased (i)
P(i) unbiased(∆Ω)
)
〈e[β∆Vumb-f(i)]〉biased (i)
(13)
where β ) 1/kBT, kB is Boltzmann’s constant, and f(i) is a freeenergy shift applied to each biasing window in order to make the frequency distribution continuous across all windows. Although more complex methods, such as the weighted histogram analysis method,25 could be employed to determine the f(i) in eq 9, here, we use a simpler approach to determine f(i) for a continuous, one-dimensional free-energy curve. Namely, the energy shifts f(i) are determined by an iterative process that minimizes the differences between two overlapping potentialof-mean-force curves in the region where both are defined. The P(i) are combined to give the full, unbiased distribution as follows (i) (i) Punbiased(∆Ω) ≡ P(i) unbiased(∆Ω) ∆ωlower e ∆Ω < ∆ωupper
(14) (i+1) where ∆ω(i) upper is set to ∆ωlower whenever two windows overlap. This unbiased frequency distribution is then normalized.
3. Application: I2 in Liquid Xe 3.1. Simulation Details. To illustrate the method, a single I2 solute molecule was simulated within a solvent of 255 Xe atoms, using periodic boundary conditions with the minimum image convention.23 Standard Lennard-Jones parameters were used to
Morse26
σ (Å) (K) De (eV) R (Å-1) r0 (Å)
I-Xe
Xe-Xe
3.94 323.7
4.10 221.6
I-I
1.5417 1.864 2.666
describe solute-solvent and solvent-solvent interactions.23 The intramolecular I2 interaction was described by a Morse potential.26 Potential parameters are given in Table 1. Liquid Xe was modeled at a density of 2.5 g/cm3 or F* ) F/σ3Xe ) 0.79 in a cubic box of length 28.16 Å. Each trajectory consisted of a classical equilibration stage of 750-4500 ps with a time step of 2.5 fs, followed by a mixed quantum-classical warm-up stage of 220 ps with a 2 fs time step, followed by a 5 ns data collection stage with a time step of 2 fs. Trajectories were generated using constant energy (NVE) molecular dynamics with velocity rescaling to maintain a temperature of ∼244 K during the first 375 or 750 ps of the initial classical warmup and in the first 20 ps of the mixed quantum-classical warmup stage; no velocity rescaling was used in the data collection stage. Trajectories were propagated in the n ) 0 ground state, and for the purpose of efficiency, only the lowest three vibrational states are included when diagonalizing the quantum Hamiltonian to solve the vibrational Schro¨dinger equation at each time step.13 Trajectories were propagated using 12 different umbrella potentials as defined in eq 7, with ∆ω(i) umb ) -20, -15, ..., 35 cm-1 and a force constant of k ) 45.56 (1/cm-1), unless otherwise noted. For each umbrella potential, two trajectories were propagated with different classical warm-up times. Twentyfour trajectories were also propagated with no umbrella potential in order to compare the statistics obtained with umbrella sampling to a comparable quantity of data generated without any biasing potential. The reported mean and variance for each variable of interest was obtained by simple block averaging. Each trajectory was divided into 5 blocks, and block averages were generated over 10 blocks obtained from 2 trajectories. For the unbiased trajectories, 12 different block averages were obtained from the 24 trajectories, and the resulting 12 averages were again blockaveraged to obtain the final averages and error bars. Error bars were calculated at a 95% confidence level using the Student t distribution.27,28 3.2. Frequency Distributions. The frequency distributions obtained with different sampling methods are shown in Figure 1. The vertical axis is plotted logarithmically to show the distribution tails more clearly. Shown in black is the composite distribution obtained from 12 umbrella sampling evenly spaced bias potentials as described in section 3.1. The smooth curve and minimal error bars indicate that sufficient sampling is achieved to resolve the distribution of line shifts ranging from -20 to 40 cm-1. For comparison, distributions obtained without umbrella sampling are displayed on the same graph. As expected, all distributions are the same within error. Without umbrella sampling (red and blue curves), only the most probable line shifts are sampled adequately in the distribution, and the error bars increase as the populations decrease; long simulation times are required to access shifts with larger magnitudes, which occur more rarely. Line shifts between -8 and 25 cm-1 are reasonably sampled within a single data set, while the range of frequencies accessed within 12 data sets is only a few cm-1 wider. With the same total data collection time, 12 umbrella potentials with overlapping windows allow access to line shifts
316 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Morales and Thompson
Figure 1. The I2 in Xe normalized frequency shift distribution obtained with different sampling methods is plotted. Results based on 12 umbrella sampling simulations (black line), 12 unbiased simulations (red line), and 1 unbiased simulation (blue line) are shown. Error bars are shown at intervals of about 1 cm-1. Results are shown for all frequencies sampled from -25 to 40 cm-1.
from -25 to 40 cm-1 and dramatically reduce the dependence of statistical error on the absolute probability. Figure 1 shows that a -20 cm-1 shift from the gas-phase frequency is over 10,000 times less likely than a 35 cm-1 shift to the blue, yet umbrella sampling allows their populations to be determined with comparable accuracy. 3.3. Mechanistic Analysis. In a previous study,3 we addressed several questions about the mechanism of vibrational frequency shifts in the I2 (and ICl) in liquid Xe system, including (1) how many solute-solvent interactions determine the line shift distribution and (2) where the solvent atoms are with the strongest effect on the line shift relative to the solute. In this section, we examine how umbrella sampling makes such calculations significantly more efficient, thereby allowing the clarification and extension of such mechanistic analyses. Beginning with the first question, we note that the explicit contributions of each solvent atom j to the I2 fundamental line shift, ∆ω ) ω - ω0, can be calculated at each time step. Specifically, these are given by ∆V(j) 01/p according to eq 5, and these are recorded at every 10th time step in the simulations and ranked in descending order 1) 2) N-1) ∆V(Rank g ∆V(Rank g ... g ∆V(Rank 01 01 01
(15)
Because the relative magnitudes of these successively ranked solute-solvent interactions are line-shift-dependent, we define a line-shift-dependent average of the “effective contributions” from the ranked solvent atoms. That is, at every 10th time step, the results are sorted by the total instantaneous line shift into histogram bins that are 0.5 cm-1 wide. Within each bin, the ensemble average is weighted so as to remove the effect of the biasing potential, V(i) umb
Average Effective Contribution for Atom Ranked k ) (k-1) 1 k (i) (Rank i) i) ∆V01 ∆V(Rank eβ∆Vumb 01 p i)1 i)1
〈 (|∑
| |∑
|) 〉
L
(i)
〈eβ∆Vumb〉L where
L - 0.25 cm-1 < ∆ω - ∆ω(0) < L + 0.25 cm-1 (16)
Figure 2. The top three effective contributions to the absolute value of the line shift, plotted as a function of the total solute-solvent N-1 i) contributions |∑i)1 ∆V(Rank |. The absolute value of this sum 01 N-1 (Rank i) |∑i)1 ∆V01 | is also shown (dashed line). The unbiased probability distribution for ∆ω - ∆ω(0) is shown (dot-dashed line) to emphasize the range of line shifts accessed. Average values and error bars for the effective contributions were obtained from 12 umbrella sampling simulations, as described in section 3.1.
The results from this analysis are shown in Figure 2, where the line shift dependence of effective contributions from the top three solvent atoms is plotted. Umbrella sampling allows these to be resolved for arbitrary line shifts; here, we focus on data for line shifts between -20 and 30 cm-1. The unbiased probability distribution of ∆ω - ∆ω(0) serves as a backdrop to emphasize the more limited range of conventional sampling. In a previous study,3 it was found that red shifts were approximated by the single strongest contribution from one solvent atom, but at least two solvent atoms are necessary to obtain a qualitatively accurate estimate of strong blue shifts from the gas-phase frequency. Poor sampling of red shifts larger than a few cm-1 obscured the possibility that strong red shifts may also be the result of more than one solvent atom contribution, that is, the error bars for the average effective contribution of the first-ranked atom were too large to distinguish it from the total sum for unbiased simulations of a total time of 24 ns. From Figure 2, it is clear that the strongest single contribution is not, in general, equal to the total solute-solvent contribution, which is indicated by the dashed line. For strong line shifts of (20 cm-1, the second solvent atom is shown to contribute about 15%, and the third solvent atom also has a nonzero effect on the magnitude of the shift. For both red shifts and blue shifts, the average effective contribution of the second solvent atom increases with the line shift. While strong line shifts in either direction cannot be explained without simultaneous interactions between I2 and two or more solvent atoms, Figure 2 shows that for red shifts of -5 to 0 cm-1, the contribution from the second solvent atom is no more than ∼1 cm-1. This is important because the unbiased distribution of line shifts at the simulation temperature, shown as a dot-dashed line in Figure 2, gives a low probability for red shifts stronger than -5 cm-1. To fully understand the role of particular mechanisms for frequency shifts and, particularly, the trends with line shift, sufficient sampling is required over a wider frequency range than that of the unbiased distribution. This example illustrates the quality of vibrational frequency data that can be obtained from umbrellasampled molecular dynamics simulations. By providing a clear prediction of the mechanisms as a function of frequency shift, such data can serve as a benchmark for theoretical models of spectral shifts.
Umbrella Sampling of Solute Vibrational Line Shifts
J. Phys. Chem. B, Vol. 112, No. 2, 2008 317
Figure 3. Spatial probability distributions for the location of the top contributing solvent atom at I2 line shifts of -2, 0, and 5 cm-1; data are collected for line shifts within (1 cm-1 of the reported value. The horizontal and vertical axes represent the distance (in Å) from the I2 centerof-mass parallel and perpendicular to the bond axis. The van der Waals contact radius and the minimum of the Lennard-Jones potential well for Xe-I2 (solid black lines) are shown for reference. Contours for red-shifting solute atom contributions to the line shift are shaded from pink to red, while blue-shifting contributions are shown in light to dark blue. The first and third columns report data obtained without umbrella sampling in 150 ps and 120 ns of simulation time, respectively. Data for the plots in the second column are from three 50 ps umbrella-sampled trajectories with an umbrella potential force constant of k ) 455.6 (1/cm-1).
The second question about the mechanism involves the locations of solvent atoms responsible for the top contributions to line shifts of a particular sign and magnitude. This is presented in the form of spatial distributions, which are easily visualized in contour plots as shown in Figures 3 and 4. In each plot, locations of contributing solvent atoms with respect to the diatomic I2 solute are compiled as a two-dimensional histogram for total line shifts within (1 cm-1 of a given value according to
Spatial Probability (xi, yj) for Atom Ranked k ) (i)
(i)
〈δ(x - xi)δ(y - yi)eβ∆Vumb〉L/〈eβ∆Vumb〉L where
L - 1.0 cm-1 < ∆ω - ∆ω(0) < L + 1.0 cm-1
(17)
In the contour plots, red-shifting contributions are indicated by shading from pink to red (corresponding to a negative probability scale), while blue-shifting contributions are indicated in shades of blue. We note that individual contributions of opposite sign may combine to give a total shift of specified sign and magnitude. For example, even when the overall line shift is 35 cm-1 toward the blue, there is still some probability that the role of the second rank atom is to shift the frequency toward the red. The level of detail in these contour plots places a premium on sampling within the selected frequency range. As a guide to the eye, the van der Waals radius and LennardJones well minimum are also depicted for the solute-solvent interaction potential with I2 at its equilibrium bond distance. The spatial distributions obtained with umbrella sampling can be examined in two contexts, (1) at total line shifts that are within the envelope of the distribution of instantaneous frequencies (see Figure 2) and (2) at line shifts that are outside of the envelope. In this way, the potential utility of the umbrella
sampling approach to obtaining mechanistic information at frequencies that are and are not sampled in an unbiased simulation can be evaluated. Naturally, the former case is one in which umbrella sampling may not be necessary. On the other hand, umbrella sampling may be able to accelerate the mechanistic analysis, particularly if fine resolution of the line shift dependence is desired. The comparison of unbiased and umbrella sampling results for line shifts within the envelope of the distribution of instantaneous frequencies is presented in Figure 3. The converged spatial distributions shown in the right-most column of Figure 3 illustrate how the spatial distribution of the top contributing solvent atom changes with the overall line shift ∆ω. For 5 cm-1 line shifts (top row), the top contributing solvent atom is situated toward one end of the I2 molecule, near the van der Waals contact radius. The atom pushes up against the repulsive wall, compressing the I2 solute along the bond axis. For red shifts of around -2 cm-1 (bottom row), the top contributing solvent atom is situated in one of two locations. Most of the time, this atom is situated in a T-shaped geometry with respect to I2 and located near the van der Waals contact radius (outlined in black). At this position, the solvent atom samples the repulsive potential wall, pushing the solute atoms apart to induce a red shift. In the other motif, the solvent atom is positioned toward one of the ends of the I2 solute but sits outside of the minimum of the Lennard-Jones potential well. It effects a red shift by pulling the solute atoms apart through attractive interactions. Not surprisingly, spatial distributions for the top contributor to line shifts of 0 cm-1 (middle row) show a combination of these three main types of solute-solvent interactions. Qualitatively, this information can be obtained from three umbrella-sampled trajectories with much shorter simulation times of 50 ps, as shown in the middle column of Figure 3. Here, a stronger force constant of 455.6 (1/cm-1) was used so
318 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Morales and Thompson
Figure 4. Spatial probability distributions for the location of the top two solvent atoms contributing to I2 line shifts of -20, -5, 10, and 35 ( 1 cm-1. Data for these plots were obtained from four 10 ns umbrella-sampled data sets, with a force constant of k ) 45.56 (1/cm-1).
that for ∼90% of each trajectory, line shifts within (1 cm-1 of the desired value were sampled. For comparison, spatial distributions from 150 ps of unbiased simulation time are shown in the left-most column of Figure 3. Distributions obtained at 0 and 5 cm-1 are similar to those obtained with umbrella sampling, but the distribution at -2 cm-1 (bottom left) is not as clear, lacking appreciable sampling of end-on attractive interactions and showing evidence of a blue-shifting, end-on repulsive interaction that is not visible in the converged histogram. The advantage of umbrella sampling is that specific line shifts are targeted to quickly obtain mechanistic information; with unbiased sampling, there is no guarantee that such information will be accessed in such a short time. Moreover, it is likely that the umbrella sampling approach may be of greater utility for systems with a broader distribution of frequencies than that for I2 in Xe. The spatial distributions for line shifts of larger magnitude cannot be reliably obtained from unbiased simulations within 120 ns. Thus, the results shown in Figure 4 represent new information that can be provided by umbrella sampling. Solvent distributions for the top two atoms contributing to larger shifts obtained from umbrella sampling simulations are shown in Figure 4. For red shifts of -5 cm-1 and stronger, both atoms have a red-shifting effect, and the T-shaped motif is strongly favored over an attractive, end-on interaction. For blue shifts of 10 cm-1 and stronger, the first atom is found in an end-on repulsive interaction with the solute, and the second atom has a ∼95% probability of being in a similar end-on configuration. For I2 to reach a 35 cm-1 blue shift, the first solvent atom must be nearly collinear with the bond axis and is located within the van der Waals radius drawn for I2 at the equilibrium geometry.
We note that within a mixed quantum-classical simulation, the I2 bond distance is allowed to fluctuate; here, the solvent is in contact with an I2 solute that is already compressed, probably by an interaction with the second solvent atom. Conversely, I2 red shifts of -20 cm-1 involve a solvent atom that penetrates the van der Waals radius at a position very close to the centroid of the I2 bond, while the second-ranked solvent atom does the same or pulls the I2 bond by an end-on attractive interaction. Such cooperative interactions involving two solvent atoms are clearly observed for strong line shifts that are rarely sampled in an unbiased simulation but are easily accessed with umbrella sampling. Table 2 lists the relative errors in the bin populations of these two-dimensional spatial histograms. These are defined over all bins j that have population Pj > 0.005, with error bars (δPj at the 95% confidence level, as follows
x∑ 1
Relative RMS Error )
M
M j)1
1
M
δP2j (18)
∑ Pj
M j)1
With conventional sampling, the statistical error based on 12 data sets (each corresponding to two 5 ns trajectories) depends strongly on the magnitude of the line shift. The RMS errors obtained by conventional sampling are less than those obtained by umbrella sampling only in the spatial histograms corresponding to smaller overall line shifts, 0 or 5 cm-1. Note, however,
Umbrella Sampling of Solute Vibrational Line Shifts
J. Phys. Chem. B, Vol. 112, No. 2, 2008 319
TABLE 2: Relative Root-Mean-Square (RMS) Errorsa in Populated Spatial Histogram Binsb solvent rank
-1 ∆ω(i) umb/cm
umbrella
no biasc
first first first first first first first first first first first first second second second second second second second second second second second second
35 30 25 20 15 10 5 0 -5 -10 -15 -20 35 30 25 20 15 10 5 0 -5 -10 -15 -20
0.07 0.07 0.07 0.07 0.06 0.07 0.07 0.06 0.09 0.08 0.10 0.10 0.06 0.06 0.06 0.06 0.05 0.06 0.06 0.03 0.07 0.08 0.08 0.08
2.02d 1.70 0.82 0.50 0.13 0.05 0.03 0.02 0.33 1.48d 2.26d 2.15 0.98 0.10 0.06 0.06 0.03 0.01 0.09 2.41d
a Root-mean-square error divided by the average population in the spatial histogram bins with populations greater than 0.005. b Total population integrated over all bins is normalized to unity. c Calculated from 12 data sets, except where otherwise noted. d Calculated from six data sets with nonzero populations.
that the umbrella sampling results for the 0 and 5 cm-1 line shifts each come completely from the two 5 ns trajectories with the biasing potential equal to that value. Thus, we see that the RMS error is only a factor of 2 larger, with effectively 1/12 of the total trajectory time. Errors are comparable for line shifts of 10 ( 1 cm-1 and more than a factor of 2 larger than the umbrella sampling errors for all other line shifts. Moreover, strong red shifts of -15 or -20 ( 1 cm-1 are not sampled at all in the unbiased trajectories, and only 6 of the 12 data sets contain red shifts of -10 ( 1 cm-1 or strong blue shifts of 35 ( 1 cm-1. With umbrella sampling, RMS errors are consistently within 10% of the average bin populations for all spatial histograms so that the solvent positions responsible for specific line shifts can be compared on an equal footing. 4. Conclusions An umbrella sampling method based on the frequency of a solute vibration has been presented. This approach does not require the calculation of additional gradients and, thus, does not substantially increase the computational cost or algorithmic complexity. Although the mixed quantum-classical approach is emphasized here, frequency umbrella sampling can, in principle, be applied within the standard perturbation theory approach, as discussed in the Appendix. To illustrate the approach, the frequency umbrella sampling has been applied to mixed quantum-classical molecular dynamics simulations of the I2 in Xe system to sample over a wide range of line shifts, as well as over a selected narrow range. With this improved sampling, spatial distributions were constructed that show how two solvent atoms act cooperatively to produce large total shifts of the I2 vibrational frequency. Additionally, because this method decouples statistical error from the magnitude of the line shift, trends in effective contributions to the line shift from successive solvent atoms were clarified for stronger shifts. Finally, spatial distributions
for the top contributing atoms were constructed over narrow ranges of frequency shifts from trajectories of just tens of picoseconds, using umbrella sampling to bias the dynamics toward the desired frequency range. Simulations like the ones presented in this work provide a view into the mechanisms by which solvent motions influence vibrational frequencies in condensed phases. With detailed line-shift-dependent data, mechanistic models can be tested in the limit of large and small instantaneous line shifts. In this context, frequency umbrella sampling is useful as a means to efficiently sample configurations associated with specific vibrational frequencies. Furthermore, it offers a clearer picture of the mechanistic motifs by providing a window into the solvent configurations that influence the vibrational frequency, for example, the cooperative influence of two solvent atoms on the I2 vibration. This ability to sample rare solvent configurations may be increasingly important for modeling vibrational spectra in complex environments, for example, that of a probe molecule in a nanoconfined solvent. For example, vibrational spectra of small molecules such as N3-,29 CH3I,30 CS2, and CHCl331 have been monitored as probes of liquid structure and dynamics in confinement. The vibrations of triatomics are amenable to mixed quantum-classical simulations, while perturbation theory can be used for larger molecules. For such heterogeneous systems, sampling over all distinct configurations is inherently more challenging, and an umbrella sampling method is a significant asset. Such spectra, and their mechanisms, are the focus of our current investigations. Acknowledgment. W.H.T. is grateful to Casey Hynes for his generous support and encouragement over many years. This work was supported by the National Science Foundation (Grant CHE-0518290). The authors thank Dr. Tolga S. Gulmen for several useful discussions. Appendix In a perturbation theory simulation, a reference Hamiltonian is chosen that is easily solvable. For convenience, we assume that this is taken to be the gas-phase Hamiltonian for the vibration. We note in passing that the solute bond is frozen at the equilibrium bond length throughout the simulation. The reference Hamiltonian hˆ (0) r is used to obtain a set of zerothorder energies and eigenvalues (0) (0) (0) hˆ (0) r φn ) En φn
(19)
Vibrational energies are then obtained as follows, to first order (0) (0) ˆ r - hˆ (0) En ) E(0) n + 〈φn |h r |φn 〉
(20)
The vibrational frequencies are then straightforwardly obtained at each time step from the energy gap as
p∆ωmn ) ∆Emn - ∆E(0) mn (0) (0) (0) ) 〈φ(0) ˆ r - hˆ (0) ˆ r - hˆ (0) m |h r |φm 〉 - 〈φn |h r |φn 〉
(21)
An umbrella potential ∆V(i) umb can be applied, as in eq 6, for a given force constant k and chosen line shift ∆ω(i) umb, resulting in the following biased classical Hamiltonian
1 2 (i) H ˜ (i) ) H0 + ∆V(i) (22) umb ) H0 + k[p(∆ωmn - ∆ωumb)] 2
320 J. Phys. Chem. B, Vol. 112, No. 2, 2008 Here, H0 is the classical Hamiltonian of the solute-solvent system with the solute frozen at its equilibrium bond distance(s). The equations of motion for the biased system are then given by equations completely analogous to eqs 9-12. In this case, however, the gradients are evaluated using the zeroth(0) order wave function φ(0) ˆ /∂e)| n , for example, (∂En/∂e) ) 〈φn |(∂V (0) φn 〉, as opposed to using the adiabatic wave function in the mixed quantum-classical approach. The standard umbrella sampling equations (eq 13) are then applied to recover unbiased properties. References and Notes (1) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698-1702. (2) Li, S.; Thompson, W. H. Chem. Phys. Lett. 2005, 405, 304-309. (3) Morales, C. M.; Thompson, W. H. J. Phys. Chem. A 2007, 111, 5422-5433. (4) Deleted in revision. (5) Fischer, S. F.; Laubereau, A. Chem. Phys. Lett. 1975, 35, 6-12. (6) Davis, P. K.; Oppenheim, I. J. Chem. Phys. 1972, 57, 505-517. (7) Oxtoby, D. Mol. Phys. 1977, 34, 987-994. (8) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187-199. (9) Li, S.; Thompson, W. H. Unpublished results. (10) Oxtoby, D. W.; Levesque, D.; Weis, J.-J. J. Chem. Phys. 1978, 68, 5528-5533. (11) Oxtoby, D. W. AdV. Chem. Phys. 1979, 40, 1-48. (12) Oxtoby, D. W. AdV. Chem. Phys. 1981, 47, 487-519. (13) Li, S.; Thompson, W. H. J. Phys. Chem. A 2003, 107, 8696-8704.
Morales and Thompson (14) Thompson, W. H. J. Chem. Phys. 2003, 118, 1059-1067. (15) Echave, J.; Clary, D. C. Chem. Phys. Lett. 1992, 190, 225-230. (16) Wei, H.; Carrington, T., Jr. J. Chem. Phys. 1992, 97, 3029-3037. (17) Colbert. D. T.; Miller, W. H. J. Chem. Phys. 1992, 96, 19821991. (18) Dickinson, A. S.; Certain, P. R. J. Chem. Phys. 1968, 49, 42094211. (19) Meyer, R. J. Chem. Phys. 1970, 52, 2053-2059. (20) Light, J. C.; Hamilton, I. P.; Lill, J. V. J. Chem. Phys. 1985, 82, 1400-1409. (21) Lanczos, C. J. Res. Natl. Bur. Stand. (U.S.) 1950, 45, 255-282. (22) Saad, Y. Numerical Methods for Large EigenValue Problems; Halstead: New York, 1992. (23) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (24) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (25) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1339-1350. (26) Herzberg, G. Molecular Spectra and Molecular Structure, 2nd ed.; Krieger Publishing Company: Malabar, FL, 1989; Vol. 1. (27) Shoemaker, D. P.; Garland, C. W.; Nibler, J. W. Experiments in Physical Chemistry, 5th ed.; McGraw-Hill: New York, 1989. (28) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, 9th ed.; Dover: New York, 1970. (29) Sando, G. M.; Dahl, K.; Zhong, Q.; Owrutsky, J. C. J. Phys. Chem. A 2005, 109, 5788-5792. (30) Czeslik, C.; Kim, Y. J.; Jonas, J. J. Raman Spectrosc. 2000, 31, 571-575. (31) Nikiel, L.; Hopkins, B.; Zerda, T. W. J. Phys. Chem. 1990, 94, 7458-7464.
