J. Phys. Chem. 1993,97, 8864-8873
8864
ARTICLES Optimal Molecular Control in the Harmonic Regime: The Methylene Halide Chemical Series and Fluorobenzene Charles D. Schwieters, Kenneth Yao, and Herschel Rabitz' Department of Chemistry, Princeton University, Princeton. New Jersey 08544 Received: February 23, 1993; In Final Form: June I , I993
In this paper, optimal control theory is applied to molecules described by harmonic force fields for the purpose of selectively exciting vibrational motion. The molecular parameters are taken from experiment, where possible. The results for different targets in the methylene halide chemical series are studied to learn how the control results track with the changing halogen constituent. It is seen that large differences in control results for different halogen constituents occur even when the target is a CH stretch. The in-plane modes of fluorobenzene are also studied to illustrate control in a larger realistic molecule with more complex connectivity. Aspects of molecule orientation relative to electric field polarization are addressed. In addition, a novel method of displaying dynamical motion is employed. I. Introduction Recently, there has been theoretical activity regarding the possibility of controlling the motion of internal degrees of freedom of a molecule by use of a laser pulse with a relatively arbitrary time dependence.1-4 An ultimate goal of the work includes obtaining laser selective chemical reactions. Previous work on this subject includes control of vibrational motion described by both quantum mechanicslJ and classical mechanic^,^ control of rotational m o t i ~ nand , ~ control of electronic degrees of freedom.4 For the special case in which (1) the potential and kinetic energy of the molecule are represented as quadratic functions of displacement and momentum, respectively, (2) the moleculefield interaction is linear in both displacement and field amplitude, and (3) the requested molecular response has a special form, there exists a closed form solution for the electric field. This solution, in addition to providing a computationally inexpensive method to obtain electric fields, allows some insight into the control process. Since the only complete potential data for many larger molecules are harmonicforce constants and nonlinear dipole data are rarer still, it is appropriate to study larger molecules using this restricted approximation to gain further intuition about the molecular control process. The goal of this study is to learn qualitative trends which might be expected to hold for more realistic potentials. Toward this goal, we examine control results for a chemical series given the same control criteria where we might expect to see similarities as well as differences due to increasingatomicweight and bonding characteristics. Other aspects we address in this study are the effects of having many modes with nonzero dipole derivatives, the effect of more complex molecular connectivity, and the role of molecular orientation upon control. In section 11, we summarize the optimal control theory formalism and introduce the methods used in the examples. In section I11 we study the optimal control results for two different targets in the methylene halide series, CH2F2, CHzC12, and CH2Br2, and in section IV we examine the control results for fluorobenzene. Finally, section V contains a brief summary.
II. Formalism and Implementation Notes In this paper, the molecular potential is taken to be quadratic in molecular displacement coordinates (bond stretches and bends).
In addition, the kinetic energy is taken to be quadratic with respect to the momenta canonically conjugate to these displacements. This latter specification is an approximation insofar as the transformation from rectilinear coordinates (where the kinetic energy has a purely quadratic dependence on momentum) to curvilinear internal coordinates is nonlinear. Under these approximations, one can write the isolated molecular Hamiltonian as
+
Ho = 5pTGp 1 p1 T F p
(11.1)
where G is a constant matrix depending on geometry with units of inverse mass and F is a matrix of force constants.6 Here q is a vector of molecular displacements and p the conjugatemomenta. If a linearly polarized electric field is present, the complete Hamiltonian (using a classical approximation for the electric field and a linear dipole approximation for the interaction) is
HT = Ho - qT&;(f)
(11.2)
where 6 is a set of three vectors [b,,b,,b,] which correspond to dipole derivatives in the three Cartesian spatial directions and ;(t) is a time-dependent electric field. Under these approximations, the equations of motion for the quantum and classical values of q(t) and p(t) are equivalent by Ehrenfest's theorem and are given by q = Gp
p
-Fq
+K(t)
(11.3) (11.4)
In the optimal control formalism, a molecular objective is specified by introducing a cost functional, J. In this paper, we take the cost functional to have the following form:
25 = c d t [ z T ( t ) W Mz ( t ) + oea2(t)] where z is the state vector,
ZT
(11.5)
= [qT,pq, and Tis the length of
0 1993 American Chemical Society
Optimal Molecular Control in the Harmonic Regime
n
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8865
f7
bond energy of each coordinate is given by (11.7)
for i = 1, ..., N where N is the number of internal degrees of freedom. To specify a molecular objective at the final time, we introduce a constraint on the state vector at the target time:
73 72
/
HTz(T) = h
w
Y
Figure 1. Template for the dihalomethane family and the definition of the local coordinates.
TABLE I: Symmetry and Internal Coordinates I, for the Methylene Halides sym species
a ro
sym coordinate
11
j
= 2.067 au.
where ii is a vector of constraint values with dimension M,the number of constraints. The M X N matrix H selects the elements of the state vector to match up with their constraints. The problem is thus to fiid an electric field producing a response in the molecule satisfying the final time constraint, eq 11.8, which simultaneouslyminimizes the cost functional, eq 11.5, given some initial state of the molecule, (z(0) = q). There exists a unique solution to this problem if the molecular degrees of freedom are controllable, which can be expressed in a closed algebraic form presented in previous work.' To visualize the resulting dynamics of a molecule, we may plot the displacement coordinates of each degree of freedom as a function of time. However, it is difficult to picture exactly how each atom is moving using this type of plot. The results of excitation of bending degrees of freedom are not obvious, nor can one easily picture energy transfer from one part of the molecule to another. To remedy this situation short of publishing movie frames, we use ORTEP software1,to create a three-dimensional plot of a molecule during some time interval with ellipsoids at the time averagepositions of each atom representing the displacement variancetensor over the time interval. The displacement variance tensor of the ith atom is given by (11.9)
TABLE Ik IR Frequencies of CHS2 vi (cm-1) 2934.4 1512.1 1116.1 527.2 3014.1 1178.3 1437.3 1091.8
(11.8)
3003.5 1456.2 717.6 283.2 3040.4 896.4 1268.2 758.1
3016.2 1366.3 592.7 171.8 3073.8 812.7 1192.6 65 1 .O
a
symspecics
j
CHzFz
CHzClz
CHzBrz
z
Ai
X
BI
Y
B2
1 2 3 4 5 6 7 8
-0.113 0.936 -0.034 0.097 -0.195 0.014 -0.025 -1.419
-0.070 -0.319 0.046 0.023 -0.069 -0.032 0.074 -0,908
0.018 -0.502 0.025 -0.025 -0.128 -0.156 0.112 -0.382
the control interval. W is a matrix weighting molecular disturbances during the entire control interval, and we is a weight on the fluence of the electric field, while e ( t ) is the amplitude of the electric field. We define the M matrix to be (11.6)
where Fd and Gd are the diagonal parts of F and G, respectively. This choice of weights for integral excitation of the molecule corresponds to a weighted sum of "bond energies", where the
where tl and t2 are the bounds of the time interval and m and n each represent the three Cartesian directions x, y , and z. ml(t) is the position of the ith atom in the direction given by m,and Miis thevalueof the position averagedover the whole timeinterval. The plots produced by this procedure thus show an average of the molecular motion through in addition to the atomic position averaged over the time interval, &i. Analogous to the control calculations, motion ofthecenter of mass and the rotational degrees of freedom are removed from these plots. Only the ORTEP plots will be shown with the illustrations, although the discussion is also based on examining bond energy and displacement plots versus time. For the control calculations shown below, we will consider just the small distribution of molecules whose orientation is fixed with respect to the laser polarization. Because our pulses are limited to 0.2 ps and the fastest time for rotation (of CHzFz at room temperature) is on the order of 5 ps, the molecules will be considered fixed in space for the duration of the pulse. While it is possible, in principle and in practice, to orient or align molecules in the gas phase, the orientation experimentsare nontrivial in themselves. The molecules not oriented at our reference orientation will undergo a different (possibly undesirable) system response, thereby reducing selectivity. However, if the observable of interest is bond energy, the average value may be quite similar to that calculated for a molecule at the reference orientation. There is considerable flexibility in optimal control theory to obtain excellent results by 'tuning" the control weights and objective terms for any particular molecule and target. This tuning is not performed in this paper, as we are interested here in the observed trends in moving through a series of similar molecules given the same form of the cost functional. II.1. MethyleneHalides. In this subsection, we present control
8866 The Journal of Physical Chemistry, Vol. 97, No. 35, 1993
Schwieters et al.
1
-0.15
a)
0.00
I
I
I
0.05
0.10
0.15
0.20
Time (ps)
Frequency (cm-
Figure 2. Optimal control of the rl C-H stretch of methylene fluoride. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps.
results for the methylene halide chemical series. Previous experimental studies have determined the harmonic force fieldg and the dipolederivativeslo for all of the methylene halide family, except CH24, for which the dipole derivatives have yet to be reported in the literature. Therefore, this study is confined to the first three members of the methylene halide family: CH2F2, CH2Cl2, and CH2Br2. Local coordinates are defined in Figure 1, and the symmetry coordinates used in defining the dipole derivatives are shown in Table I. The molecular geometry parameters are given in ref 9. For convenience, the fundamental frequencies and dipole derivatives are given for the methylene halide family in Tables I1 and 111, respectively. Although symmetry coordinates are most convenient for succinctly specifying the molecular parameters, local bond coordinatesare more useful for the control calculationsperformed here. Both coordinate systems describe the same dynamics, but the cost functionals differ if only diagonal weighting is used in the W matrix. As seen in previous work? it was found here that the best results were obtained by using the coordinate system in which the target is most naturally expressed. The set of internal coordinates used in the control calculations shown below is listed in Table I, next to the symmetry coordinates. Note that there are only four coordinates containing bending degrees of freedom in the table, although six are shown in Figure 1. One bending motion belongs to the A2 irreducible representation, and, thus, its excitation is forbidden by a linear dipole interaction. The other unaccounted for bending motion is given by a redundancy condition. For the examples shown below, the electric field is taken as polarized at 54.74O to each of the positive x-, y-, and z-axes, with the orientation of the molecule being that shown in Figure 1: the symmetry axis is along the z-axis, while the hydrogen atoms lie in the x-z plane. This choice of polarization breaks the molecular symmetry and allows excitation of all IR active modes.
General trends in the series can be seen in the normal mode frequencies and in the dipole derivatives. All of the frequencies associated with the halogen degrees of freedom (vj, vq, V6, UT, and us) decrease down the series due to the increase of halogen mass and a decrease in bond strengths, while the frequenciesassociated with the hydrogen degrees of freedom do not change appreciably. The dipole derivatives generally decrease in magnitude down the series, due to the fact that the charge separation decreases because of the decreasing electronegativity of the halogen atom. Because all IR active modes have nonzero intensity for these molecules and there is no true degeneracy, any target state can be attained due to controllability arguments.* II.1.1. StretchingOne CH Bond. The first target studied was that of stretching one CH bond (rl), while minimizing the disturbance to the remainder of the molecule. The constraint at the target time was rl( 2') = 1 au, which correspondsto stretching one C H bond from its equilibrium value of 2.069, 2.048, and 2.048 for CH2F2, CH2C12, and CH2Br2, respectively. The weighting parameters were ~ = [ O1 1 1 1 1 1 11
4 = [ 0 0 0 0 0 0 0 01 we =
1
T = 0.2 ps where wz and 4 are the diagonal elements of the W matrix associated with the displacement and momentum vectors, respectively. The off-diagonal elements are taken as zero here and throughout the remainder of the paper. This choice for W corresponds to weighting displacements of all the degrees of
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8867
Optimal Molecular Control in the Harmonic Regime I
I
0 10 0 05 3
2 9 Y
L
w
000
-005
-0 10
-0 15
Time (ps)
Frequency (cm- ')
Figure 3. Optimal control of the rl C-H stretch of methylene chloride. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps. 0.3
I
I
0.2 0.1
-,
0.0
-3
-0.1
u
-0.2
2
a rn
-0.3L
6
I
0 00
0,
I1 u
I 0.05
I
I
0.10
0.15
0.20
Time (ps)
0
1000
2000
3000
4000
5000
Frequency (cm-')
Figure 4. Optimal control of the rl C-H stretch of methylene bromide. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps.
freedom except the target coordinate, rl. In addition, all calculationswere made using initial conditionsof zero (i.e.,z( t=O)
= 0). The results are shown in Figures 2, 3, and 4 for CH2F2, CH2C12, and CH2Br2, respectively.
8868 The Journal of Physical Chemistry, VoL 97, No. 35, 1993 O
Schwieters et al.
3t %,
3 4
0
a m d
a,
$
a I
0 00
I
I
I
0 05
0 10
0 15
V
Y
I
0 20
Time ( p s )
Frequency ( c m - I )
Figure 5. Optimal control of the r2 C-H stretch of methylene bromide. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps.
Time (ps)
Frequency (cm-')
Figure 6. Optimal control of the R I C-F stretch of methylene fluoride. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps.
Only the rl (ZI)bond is appreciably stretched in CHzF2, while the rest of the molecule is left undisturbed, for the most part. If
the bond energy is the quantity of interest, the results are even better. The beginning of beating between the symmetric and
Optimal Molecular Control in the Harmonic Regime
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8869
L
Time (ps)
Freqency (cm-')
Figure 7. Optimal control of the RI C-Cl stretch of methylene chloride. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 p.
asymmetric CH stretch frequencies is seen in the electric field plotted in Figure 2a, which is corroborated by the wide frequency distribution around 3000 cm-1 in the power spectrum in Figure 2b. Figure 2c contains an ORTEP plot depicting the molecular motion of CH2F2 averaged over the interval 0.1-0.2 ps. The center of each ellipsoid denotes the mean position of the correspondingatom, while the semimajor axis shows the direction of greatest motion. Here the plot shows that only the H I atom undergoes appreciable motion and that most of this motion is in the direction of the rl bond axis. Thus, we see that the target is met quite well in this example. The situation is quite similar for attempting to stretch the rl bond in CH2C12, although the r2 bond and 1, and 1 8 become slightly moreexcited. Moreover, the electric field has a somewhat larger peak amplitude and fluence, as seen in Figure 3a. This latter behavior is due to the decreased dipole derivatives of the CH stretches relative to those of CH2F2. The ORTEP plot in Figure 3c shows slightly more excitation of the H2 atom but is otherwise quite similar to that in Figure 2c. While the C-H portion of these molecules is rather similar, with nearly identical fundamental frequencies and similar dipole derivatives, the optimal electric fields in Figures 2a and 3a are rather different. A simple sinusoidal field at one of the C-H stretching frequencies would excite both C-H bonds because of the large coupling between the two. Thus, the optimal control formalism is able to intelligently search out a field which causes the desired molecular response. The results for attempting to stretch the rl CH stretch of CH2Br2 are shown in Figure 4 and are qualitatively different from those of the first two members of this chemical series. The results show large excitation of the r2, Rl, 17, and 1 8 degrees of freedom, with the stretch and energy of r2 actually exceeding that of rl during most of the control interval, including at the
target time. The ORTEP plot in Figure 4c shows that all atoms undergo large amplitude motion and that the bending motion of H1is much more pronounced than the desired stretching motion. Parts a and b of Figure 4 show that while the largest frequency components of the electric field are at approximately 3000 cm-1, corresponding to the CH stretches, there are significant contributions at very low frequencies correspondingto C-Br stretching and bending. The laser polarization in this set of calculations is more along rl than r2, so one would expect rl would be easier to pump. However, the r2 dipole derivative for CH2Br2 is actually larger for this molecule/laser orientation. From Table 111, it is seen that the sign of c3pz/aSl in CH2Br2 is opposite that in the other two members of the series in this study. This is evidence either of changing bonding character or of an error in the dipole derivative determination. To attempt to reduce the r2 stretch in this case, the field pumps at many different frequencies and tries to manipulate different interfering pathways to the target. This leads to excitation of the entire molecule. Figure 5 shows the results for choosing a target of stretching r2 instead of rl in CH2Br2 using the same laser polarization. The displacement of r2 shows much more specificity. While the ORTEP plot in Figure 5c shows behavior more similar to the rest of the methylene halide results for stretching r1, there certainly is somewhat more excitation of the other modes. The electric field in Figure Sa is also more similar to those in Figures 2a and 3a, with there being only one large frequency component at around 3000 cm-1. This result shows the important effect of molecular orientation relative to the laser polarization upon the quality of the control behavior. 11.1.2. Target of Stretching One CX Bond. As further illustration of trends in the chemical series for control in the harmonic regime, we chose to attempt to stretch a single C-X bond(RI)fromitsvalueof2.612,3.360,and3.657auforCH~F2, CH2C12,and CH2Br2, respectively. To this end, we constrained
8870 The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 I
0
Schwieters et al.
I
I
I
I
I
os1
-
60 -
-
30 -
-
b) 0 0
1000
Time (ps)
2000
I
I
3000
4000
5000
Frequency (cm- ’)
Figure 8. Optimal control of the R1 C-Br stretch of methylene bromide. (a) The optimal field. (b) The power spectrum of the optimal field. (c)
ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps.
TABLE IW In-Plane Internal Coordmates for C.&F no. definition description 2-6 7-12 13 14-18 19 20 21
r2, ...,r6
R I ,...,R6 81 = 1 / d 2 ( 4 1 - 4’1) 82,
.**,86
+
+
l / d 6 ( a 1 - a2 a3 - a4 as - a6) 1/d12(2a1 - a2 - a3 2a4 - as - (Ye) 1/d2(az-a3 a5-(Y6)
+
+
TABLE V In-Plane Fundamental Frequencies and Dipole Derivatives of C&F
C-H stretches C-C stretches C-F bend C-H bends ring deformation ring deformation ring deformation
R1( 7 ) = 1.O au and used the following weights in the calculation of the optimal field:
1 2 3 4 5 6 7 8 9 10 11
3109 3090 3064 1617 1501 1260 1161 1015 997 807 511
12 13 14 15 16 17 18 19 20 21
3107 3076 1607 1456 1313 1271 1163 1067 617 387
1 2
Y Y
3
Y
4 7
Y Y
8
Y
9
Y
X X
X X X
wj=[l 1 0 1 1 1 1 13
d=[O0
0 0 0 0 0 01 we=
1
T = 0.2 ps The results for CH2F2, CH2C12, and CHzBr2 are shown in Figures 6-8. The electric field is polarized almost perpendicular to the R1 bond axis, so it would be expected that this target might be difficult to achieve. Indeed, the results for CH2Fz in Figure 6 show very poor selectivity for this target. R2 is actually excited quite a bit more than R1 at the target time t = T. The ORTEP plot in Figure 6c averaged over 0.1-0.2 ps shows large lateral motion of the F2 atom due to large amplitude motion of the composite bend Is near the final time. The electric field is large in magnitude with a very broad power spectrum spanning the range from dc to about
Bz
x
13 14
Y
15
Y
16 19 20 21
x Y Y x
X X
0.637 -0.050 -0,090 0.037 -0,065 -0.079 -0.217 -0.03 1 0.098 -0,067 0.07 1 0.196 -2.849 -0.371 0.227 0.350 0.155 0.396 0.012 0.135 -0.510
2500 cm-I. This occurs in an attempt to pump the Rl bond with a very small dipole derivative and large coupling to other modes. In another calculation not shown here, we chose a target of R2 using the same polarization. As expected from the behavior in Figure 6, the results for attempting to selectively excite the R2 stretch are much more successful. It shows again that the moleculelaser polarization is crucial to the quality of the control results, and, considering the possibility of preparing the molecular orientation, perhaps a separate previbrational control procedure should be carried out where the polarization of the laser is varied to minimize the cost functional.
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8871
Optimal Molecular Control in the Harmonic Regime
0
F@m 9. Fluorohcnzene and the definition of its local coordinates.
Figure 7 shows that the selectivity is somewhat better for stretching R I in CHZClz, despite the choice made for orientation. .~---
O ;?
I
0.00
0.05
0. I O
Time ( p s )
0.15
0.20
Decent energy localization is found in R1, but relatively large displacements of Rz and especially ISoccur. The power spectrum of the electric field in Figure 7b shows a much narrower frequency spread than in the previous case for CHzFz, while the peak magnitude of the field itself is down by a factor of more than 3. Although the ellipsoid corresponding to C11 in Figure 7c is not well aligned along the R1 bond, the ellipsoids corresponding to the three lighter atoms have their principal axes primarily aligned along theRl bond, denoting excitation of this bond. The ORTEP plot shows large displacement of the carbon and hydrogen atoms because all motion is shown relative to a stationary center of mass. Thus, because the carbon and hydrogen atoms are light in comparison with the two chlorine atoms, most ORTEP motion is shown occurring in these three lightest atoms. Figure 8 shows the results for CHZBrZ, and it is seen that the optimal field in this case excites the RZbond even less relative to the desired R I bond (the ellipsoids of C, HI, and Hz are better aligned along Rl), although the displacement of Is is still quite large. The trend from CHzFz to CHzBrZ of reduced excitation of the Rz stretch while the R1 stretch is the target is mostly due to the fact that, for thisorientation, the ratioofthedipolederivative of R1 to that of Rz for this polarization increases from 0.2 to about 7, as seen in Table 111, thus making it increasingly easy to excite the R1 bond. One reason for the improved results down the chemical series is due to the increasing X-C-X angle from 108.32O in CHzF2 to 112.9' in CHzBrz. This change causes the laser polarization to be more along the C-XI bond axis and thus leads to the larger value of the dipole derivative. Even though the motion is well localized in the displacement, it is difficult to visualize the dynamics in the ORTEP plot (Figure 8c) because the bromine atoms are much heavier than the rest of the molecule.
IO
-
I
05
E
-
0 0,
a
v)
PI
b
0
M
2
I
4
00
n
I
80.-
z
B
i 3 Y d
I
-05
-
40-
a 0
-1 0
Time (ps)
I
Az
I
b)
Frequency (cm-')
-
-0 06 a)
-0 09
-
I
I
I
-
Figure 11. Optimal control of the R3 C=C stretch of fluorobenzene. (a) The optimal field. (b) The power spectrum of the optimal field. (c) ORTEP plot of the dynamics averaged from 0.1 to 0.2 ps. (d) Displacement of three internal coordinates: the solid, dotted, and dashed curves represent Ra, R2, and R6,respectively.
In center of mass coordinates, the bromines remain relatively still while the other three atoms undergo much motion. The large excitation of ISis common to each member of the series for the current target and is reflected in the ORTEP plots by the fact that the ellipsoid corresponding to XIis not oriented along the CX bond axis. While the dipole derivative of this mode is relatively small for each member of the series, it is strongly coupled to the C-F asymmetric stretch. 11.2. Optimal Control of Cd-Ifl. In this section, we study optimal control of a larger molecule, fluorobenzene, where only the in-plane modes are considered. It was imagined that more complicated energy pathways might be possible in this molecule because it has more degrees of freedom than the methylene halides. While a harmonic force field has been reported in the literature,12there has been no complete study published of the dipole derivatives. Therefore, we calculated finite difference dipole derivatives for each degree of freedom using the STO-3G basis set in Gaussian-86 and the geometry from ref 12. While these dipole derivative calculations are not expected to give quantitative accuracy, they should be chemically meaningful. The symmetry coordinates, fundamental frequencies,and dipole derivatives for each unique coordinate are shown in Tables IV and V, and the local coordinates are defined in Figure 9. For the studies of CsHsF, we took the laser to be polarized in the x-y plane at 4 5 O to the positive x - and y-axes. 11.2.1. Target: Excite One C-H Bond. The first target was to stretch one C-H bond, rz, while disturbing the remainder of the molecule as little as possible. This bond is almost directly along the axis of the laser field's polarization, so one might expect that it would be relatively easy to excite. The control parameters used in this case were
w j = [ O 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 01 $= [ l o 10 10 10 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 we = 1
T = 0.2 ps Here the momenta were weighted instead of the displacements to prevent large spikes in the energy delivered to the molecule immediately before the target time. The momenta of the other CH stretches were weighted more strongly because they are more tightly coupled with r2 than the other modes. In addition, the weighting schemein the cost functional was modified by removing the M matrix from the integral disturbance term for these studies. The inclusion of the M matrix weights each degree of freedom by its energy contribution, thus leading to less weighting on the softer modes before multiplication by the W matrix. This is undesirable if one wishes soft modes to undergo small displacements. Figure 10 shows the optimal control results for this case of stretching one C-H bond. The r2 bond (12) is indeed the most highly excited bond in the system. It is noteworthy that 1-5 ( I s ) , the C-H bond opposite r2, has a relatively low excitation even though it is also approximatelyaligned along the field polarization (however, its dipole derivative is slightly less than that of rz, as seen in Table V). All of the C-H stretches are strongly coupled
Optimal Molecular Control in the Harmonic Regime because of their similar frequencies, so it is gratifying that the optimal control formalism found an electric field which does not cause large excitation of all of the C-H bonds. The ORTEP plot in Figure 1Oc shows displacement ellipsoids averaged over 0.10.2 ps. It is seen that the excitation of the angular degrees of freedom of three hydrogens besides Ha results in relatively large motion in Cartesian space, but it is also clear that the rz bond has the most motion in the stretching degree of freedom. Parts a and b of Figure 10 show the optimal field and its power spectrum, respectively. The electric field pumps primarily at about 3000 cm-l but also has a significant dc component. 11.2.2. Target: Excite One C-C Bond. The final example was to selectively excite the R3 C-C bond in the six-membered ring. This target was expected to be relatively difficult because the ring stretches are closely coupled with many other modes of the system. The parameters used in this control calculation were
wX=[O 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 01 ~ = [ 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 ] we = 1
T = 0.2 ps The results of this calculation are shown in Figure 11. As expected, many coordinates have large displacement over much of the control interval. However, relatively good energy localization in R3 (I,) occurs at the target time (not shown here). The electric field and its power spectrum are shown in parts a and b of Figure 11. The electric field has a quite complicated shape with many frequency components. It is seen that the peak frequency component of the optimal field is at dc, with the next largest peak at 130&1400 cm-1, where one would expect to pump to excite a C-C stretch. Surprisingly,there is also a nonnegligible frequency component at the C-H stretchingfrequency,reflecting that the whole molecule becomes excited in this example. The ORTEP plot of displacement averaged from 0.1 to 0.2 ps in Figure 1ICshows the large amplitude motion of the whole molecule. The large excitation of the C-H bends is seen in the very large ellipsoids of the hydrogen atoms. The R3bond does not appear very excited in this picture, but notice that both the C3 and Cq ellipsoids are oriented more or less along the R3 bond axis, illustrating the stretching nature of their motion. The large dc component of the field also causes the average position of the fluorine atom to be bent quite a distance from its equilibrium position. Figure 1Id shows the displacement of the three CC bonds with the largest displacement,Rz,R3, and Rg. Here it can be seen that the stretch of the R3 bond is about 20% greater than that of the bond with the next largest displacement at the target time. The reason for the poor selectivity in this example is that both the C-F stretch and the C-H bends lie in the same frequency region as the C - C ring stretches. In addition, these other modes have larger dipole derivatives, and all of them are tightly coupled, so selective excitation of the target bond is indeed a very difficult goal here. We could constrainall the coordinatesto an acceptable state at the final time; however, this would cause the results to be worse before this time.
III. Summary In this paper, we examined trends along a chemical family of molecules and found both similarities for the members of the family and differences expressed as trends while traversing the halogen column of the periodic table for the methylene halides. Differences were caused by changes in geometry and the nature
The Journal of Physical Chemistry, Vol. 97, No. 35, 1993 8873
of the C-X bonding. In addition, we examined a larger molecule, fluorobenzene,and saw that more complicated field shapes might arise when a large number of modes were coupled to the target degree of freedom. It was seen in this study that molecular orientation plays an important role in the quality of the results, and this is readily understood intuitively. If it is possible to control the orientation of the molecule with respect to the polarization of the electric field, it may be desirable to choose this orientation so that it optimizes the cost functionalfor a particular target. Recent work has explored the possibility of rotational5 and orientational13 control. This optimization step would be a nonlinear procedure but would be fairly simple, with only three degrees of freedom. It would only be necessary to maintain orientationof the molecule in space during the vibrational control interval. However, the effect of rotation onvibrationalcontrol resultsfor control intervals similar to those chosen in this paper was seen to be negligible in a previous study.14 Another consideration is the choice of the weights in the cost functional. These weights are crucial to obtaining acceptable results, and it is not necessarily obvious how to choose their values u priori. One choice is to weight all modes proportional to their bond energy, as was the case in the methylene halide study, but this can lead to very small weighting for modes with small force constants, such as bending modes, and thus large excitation in the optimal control result. For this reason, the weights were chosen independent of bond energy in the fluorobenzene example. We have introduced a method of visualizing dynamic motion using still frames by utilizing existing crystallographic thermal ellipsoid plotting software. This was seen to give a different view of the molecular dynamicsthan is obtained by plotting molecular coordinates versus time. Because of the nature of these plots, however, motion of degrees of freedom including heavy atoms in small molecules or large-scale collectivemotion is not particularly wellvisualized. However, it was seen to be quite useful for motion of atoms which contain only a fraction of the total molecular mass. Acknowledgment. We acknowledge support from the Army Research Office and the Office of Naval Research. References and Notes (1) Shi, S.; Rabitz, H.J . Chem. Phys. 1990, 92, 364. Shi, S.; Rabitz,
H. Compur. Phys. Commun. 1991,63,71. (2) Tannor, D. J.; Rice, S. A. J. Chem. Phys. 1985,83,5013. Rice, S. A.;Tannor, D. J.; Kosloff, R. J. Chcm.Soc., Faraday Trans. 2 1986,82,2423. Kosloff, R.; Rice,S.A.; Gaspard, P.; Tersigni, S; Tannor, D. J. Chcm. Phys. 1989,139,201. Tersigni, S.H.;Gaspard, P. and Rice, S. A. J. Chcm. Phys. 1990,93, 1670. (3) Schwieters, C. D.; Rabitz, H. Phys. Rev. A 1991, 44, 5224. (4) Gross, P.; Neuhauser, D.; Rabitz, H.J. Chcm. Phys. 1992,96,2834. ( 5 ) Judson, R.; Lehmann, K. K.; Warren, W. S.; Rabitz, H. A. J . Mol. Srrucr. 1990.223.425. Shen. L.: Shi. S.: Lin. C.: Littman. M.: Rabitz. H.: Heritage, J. P.; Wiener, A. M:, unpublished resulk. Shen,L.; Rabitz, 8. Phys. Chem. 1991, 95, 1047. (6) Wilson, E. B., Jr.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: New York, 1980. (7) Beumee, J. G. B.; Rabitz, H. J. Marh. Phys. 1990, 31, 1253. Shi, S.; Rabitz, H. J. Chcm. Phys. 1990, 92, 2927. (8) Schwieters, C. D.; Beumee, J. G. B.; Rabitz, H.J. Opt. Soc. Am. B 1990, 7, 1736. (9) Ford, T. A.; A r m M., R.; Robinson, E. A. S. Afr. J . Chcm. 1977, 30, 95. (10) Ford, T. A.; A r m M., R.; Robinson, E. A. Spectrochim. Acta 1978, 34A, 17. (1 1) Johnson, C. K. OR TEP-II: A FORTRAN Thermal-Ellipsoid Plor Programfor Crysral Srrucrurc Illustrations; NationalTechnical Information Service, ORNL-5138, 1976. (12) Fogarasai, G.; C h a r , A. G. Spectrochfm.Acfa 1988, I I A , 1067. (13) Parker, D. H.; Bernstein, R. B. Annu. Reu. Phys. Chem. 1989,40, 561. (14) Amstrup, B.; Carlson, R. J.; Matro, A,; Rice, S. A. J. Phys. Chem. 1991, 95, 8025.
