Maximum entropy imaging and quantum molecular timescale generalized langevin equation theory

Chemical Physics ELSEVIER

Chemical Physics 211 (1996) 91- 100

Maximum entropy imaging and quantum molecular timescale generalized Langevin equation theory H. Keith McDowell a, A.M. Clogston a Department of Chemistry and Biochemistry. Universi O, of Texas at Arlington. Box 19065. Arlington. TX 76019-0065. USA b AT& T Bell Laboratories. Murray Hill, NJ 07974-2070, USA

Received 23 October 1995; in final form 20 May 1996

Abstract

A new maximum entropy (MaxEnt) imaging procedure is introduced to treat the inverse problem in quantum molecular timescale generalized Langevin equation (MTGLE) theory. In particular MaxEnt imaging is used to construct the MTGLE bath spectral density from knowledge of its moments.

1. Introduction

Time correlation functions (TCFs) play an essential role in our understanding of the dynamics of large-scale atomic and molecular systems [1,2]. In particular one can obtain information about relaxation times and the location and shape of spectral peaks from TCFs [3-5]. Recently, we introduced a quantum version of molecular timescale generalized Langevin equation (MTGLE) theory designed as a computational tool to calculate such properties [3-5]. The method was shown to be efficacious in this regard in application to a number of systems. MTGLE theory is based on the notion that information about the TCF of interest can be obtained from knowledge of its time derivatives at time zero. In particular and following previous notation [3-5], if we denote the desired TCF as 2 ( t ) where the dot

i Retired.

indicates time differentiation, then the time derivatives are given by o-0~n)= ( - 1 ) " [ d 2n )(( t ) / d t 2~] ,=0"

(1)

We assume that the Hamiltonian describing the system is stationary and thus )~(t) is an even function of time [6]. We also assume that )~(t) is normalized such that ,{'(0)= 1 = o-0(°). We define a spectral density for 2 0 ) given by 2 P(o)) = 7 £

d t ) ~ ( t ) cos o)t.

(2)

It can be easily shown that [3] ~c

or0(")= f0

do) o)2]O(o));

(3)

that is, the time derivatives in Eq. (1) are equivalent to the moments of the spectral density in Eq. (2). The MTGLE procedure is accomplished by first taking the Laplace transform of /((t) and defining a bath TCF O~(t) such that [3] ~ ( z ) = [z 2 + w2 - w 4 , 0 , ( z ) ] - ' ,

0301-0104/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0301-0104(96)00224-8

(4)

H.K. McDowell. A.M. Clogvton / Ct en cal Pt y~'c~- 211 (1996) 91-100

92

where the Laplace transform is designated by a caret and the transform variable is z. From Eq. (4), the parameters to2co and o)4c I are easily shown to be related to the time derivatives in Eq. (1) such that :

(5)

and

to/, =

O-0( 2 ) - (O'0(')) 2.

(6)

Using Eq. (4), it can be shown that the time derivatives of Ol(t) given by o~ ~) = ( - 1)~[d2~{)1(r)/d,2~],=0

(7)

can be computed from the relationship n-2

O-(n) = toe,20.0(n

o_lJ)o_;n-2-j),

l)_}_ tool4 £

(8)

j=0

with o-~°) = 1 [3]. One also defines the bath spectral density as 2

Pl(to) = ' T f 0 dt (},(t) cos tot.

(9)

In our MTGLE approach we expand Or(t) in a set of specially designed GLE functions Gk(to c, %; t) such that [3] M+I 0 1 ( t ) : E OlkOk(toe' toe'; !7)

The process of using the moments or('') to obtain spectral densities such as p~(to) is an example of the classic inversion problem [8] and was explored in early work by Gordon [9,10]. One approach which has emerged to deal with this numerical problem is to use the maximum entropy (MaxEnt) principle [1,8]. Numerous variations of the moment MaxEnt strategy have been pursued [1,8,11 22] and arguments advanced in support of information entropy as the method of choice [23,24]. It is the purpose of this paper to examine a variation of the Gull-Daniell MaxEnt imaging approach [11] as a means of computing the bath spectral density Of(Oo) from knowledge of the moments. The need for this approach arose from our experience with other forms of MaxEnt imaging in the context of using moment information and the effect of the linear dependence among the moments on nonlinear search algorithms to determine Pl(to)We begin in Section 2 by presenting our MaxEnt approach. In Section 3 we apply the procedure to a known Bessel function TCF. Concluding remarks are made in Section 4.

2. MaxEnt imaging

(10)

To carry out a MaxEnt imaging procedure in the context of MTGLE, we define an information functional I given by

(11)

I= - -~,f' dto p,(to) lneg (`to~----~ ILL

k=l

and M+I

~,(z)

=

o~d,(to,,, to~; Z),

PI(~)

a

k=l

where toe and to, are two parameters which control the behavior of the GLE functions and M is the number of known bath moments o-~") [7]. The coefficients % are determined by requiring the expansion in Eq. (10) to have the correct moments up to M. Once the expansion in Eq. (10) is determined, the bath spectral density PL(to) can be computed from the relation Pl(to) = - ( 2 t o / T r ) I m

~'~

0,(ito).

(12)

A full description of this MTGLE process is presented elsewhere [3,7]. For the purposes of this paper we focus our attention on the computation of p~(to). Berne and Harp suggested a similar focus based on moments of the bath spectral density in 1970 [1].

- /3f¢% dto p , ( t o ) In g ( t o ) , roa

(13)

where the first term is a cross entropy [25,26], the parameter L is given by M

L =

£ , . y ( m , ( ~ , n ) _ _ O. 1( m ))2. m= 0

(14)

the parameters /x and /3 are regularization parameters 2, g(to) is a specified distribution which can be a prior distribution [15] or default model [28] or, as

2 For a recent discussion on regularization parameters, see Ref. [27].

H.K. McDowell, A.M, Clogston / Chemical Physics 211 (1996) 91 I00

in the present case, a computational aid to rapid convergence of the maximization process, the factor e is included for normalization purposes, cr~ ") are the known bath moments, and

do, o,2>,(o,).

(15)

co a

The form for the least squares function L in Eq. (14) allows for 3/(") to be chosen to keep large moments from numerically overcoming small moments. In the present work we choose y(,,l = (crtIm))-2. The integration endpoints w a and w b are assumed to be chosen so as to encompass the range required for p~(w). The information functional I is effectively a constrained cross entropy. The conventional form for the information functional in the MaxEnt literature is a S - L with the generalized Shannon-Jaynes entropy S given by S=-

dco p 1 ( w ) - m(~o) a

+p.(w)

In m ~ -

!] '

(16)

where m(a~) is a prior spectral distribution or default model, a is a regularization parameter, and the parameters ,/(,~/ in Eq. (14) are chosen to be y(") = (2cry) J with o'~ the root-mean-square error in the data [29]. In our case the data are the known moments. An excellent and concise presentation of the conventional approach to MaxEnt imaging with a Bayesian foundation is given by Gubernatis, Jarrell, Silver, and Sivia [29]. A recent application has been made by Gallicchio and Berne [30]. The functional form we define in Eq. (13) reduces to Eq. (16) for ,8 = 0 and g ( w ) = re(w) aside from a simple shift of the entropy by the normalization constant for m(~o) and scaling by /x-~ with o~=t x - I . For , 8 = I the distribution g(o)) is removed and Eq. (13) reduces to Eq. (16) with a flat default model. We have chosen the form in Eq. (13) to avoid confusion with conventional MaxEnt approaches and to achieve several purposes. First, we want to make use of the fact that MTGLE moments are often known to a high degree of accuracy. Second, we want to overcome convergence problems associated with the nonlinear search for p~(o)), which typically arises from linear

93

dependence of the moments, by using the distribution g ( w ) as a guiding distribution. Third, we want to remove any residual effect of the guiding distribution g(eo). Finally, we want to achieve an image of p~(a)) consistent with a MaxEnt approach. We elaborate on these points below. In the conventional Gull-Daniell MaxEnt process one searches for the choice of /, (or c~) which produces a spectral distribution ioj(w) which brings L to a specified value consistent with the errors in known data. For our case with L defined as in Eq. (14) and with our moments accurately known, the purpose of the process is rather different. We propose to set L at a small value and attempt to find P t(~°) such that the moments of the image are within a specified tolerance of the known values. There are many imaging schemes which can accomplish this task including the MTGLE process described in Section 1. The advantage of using Eq. (13) to obtain the image is that it satisfies a least biasing algorithm in the sense of Levine [18] as well as fitting the moments. Ultimately, the MaxEnt scaffolding used to set up the imaging process as defined by Eq. (13) and the nonlinear search process described below is not important if the images produced are routinely of a satisfactory quality. The MaxEnt strategy is to maximize I in Eq. (13) with respect to the function p~(w). Taking the functional derivative of I with respect to p~(w) and setting it equal to zero, we find that p,(,,,) = g(o,) '-e

M e~'°~- o-("~ 0, 2" ) ×exp

- 2/, Y'.

~T~

O..(m)

.

m=O

(17) In this form we see explicitly that ,8 = 1 removes the dependence on g(~o). Although the expression in Eq. (17) represents a relationship obeyed by the exact solution to the extremum problem, we have not found it to be a useful expression for direct computation since the computed moments ~ m ) in the exponential depend on pL(co) through Eq. (15). On the other hand, for a given choice of p. and ,8, Eq. (17) can used to construct a dual space numerical approach [31-33]. We will not explore this avenue at this time.

ILK. McDowell, A.,44. Clogston / Chemical Physics 211 (1996) 91 I O0

94

To obtain numerical solutions to Eq. (13). we proceed by using quadrature methods to evaluate integrals and a regularization procedure to find an appropriate value for the parameter/z. The parameter /3 is set to zero or unity in the final analysis depending on the use of g(02). In the present work we choose g ( w ) as a guiding distribution to be removed by setting /3 = 1 at the end. We note that Gubernatis et al. [29] have developed a Bayesian procedure for direct evaluation of o~ (or /x) in the standard G u l l Daniell procedure. Herein, we are interested in presenting an annealing walk on /z and 13 which can be thought of as part of the nonlinear numerical procedure to obtain Or(02). We pursue this route to avoid numerical convergence problems associated with the nonlinear search for p~(02). Defining f(02) = p~(02)/g(02), we obtain ~.v f,, I = - Y'~ G0f,, In -/xL-/3 n=l

e

x ~ r,,0f~ In g,,, n=l

(18) with N

a-~'"' = E r,,,,,J;,

£2~ 2 =

f0

do)

p,(02) - -0) 2

/*,~w = #o,~,( 1 +

in Eq. (14) for L and rnm = Wn gn 022,, n

(20)

where % is the quadrature weight, 02,, is the quadrature point, g , = g(02,), and f,, =.f(02,). In the present work we use an extended Simpson's rule to evaluate the integrals [34]. For this case we assign lower and upper frequency cutoffs to the integrals and evenly space the quadrature points with the lower and upper frequencies being the first and last points 02t and cox. If the spacing between points is denoted by h, then w~ = w N = h/3. For w,, where n is even. we have %, = 4h/3. For n odd, w,, = 2h/3. In practice we experiment with the location of the cutoff frequencies to ensure that the final results are insensitive to their chosen values. As described above, we use the distribution g(02) as a guiding distribution to enhance the nonlinear search. In the present case we use the distribution [35] ,/2 02:e ,,o-"

(21)

(22)

is finite. For diffusional problems a different choice for g(02) may be required. The Gaussian exponential decay reflects the short-time Gaussian behavior of the time correlation function 0L(t) [1,36]. In principle a MaxEnt procedure based on Eq. (13) should be numerically independent of the distribution g(02) at /3 = 1 when it is interpreted as a guiding distribution. We test this property herein by using a very accurate choice for g(02). The procedure for using the digitized version of MaxEnt imaging described by Eqs. (18) and (19) is to first set /, = / 3 = 0. Thus, we have Or(02)= g(02) as a starting point. We next increase the p a r a m e t e r / , with /3 = 0 using a jump algorithm such that

(19)

n--l

g(02) = 4 a ( a / r r )

where the parameter a is detemfined by minimizing L in Eq. (14) with Pt(02) = g(02) in Eq. (15) and the endpoints are taken to be ¢o~ = 0 and 02b = ~:- This choice for g(02) is based on the assumption that the spectral density p~(02) has a single peak and takes off from 02 = 0 with a polynomial dependence of at least 02e. This ensures that the M T G L E adiabatic frequency ~(21 for the bath given by [3]

6R),

(23)

where R is a random number between 0 and 1 and 6 is a jump magnitude. For the new choice of /, we carry out nonlinear maximization for fn in Eq. (18) using the P o l a k - R i b i e r e variant [37] of the F l e t c h e r - R e e v e s method with a convergence tolerance on the functional 1 of " f t o l " . The parameter ftol. which is defined in Eq. (27) below, is normally in the range 10 -s to 10 - ~ . The standard computer package [37] for this method has been modified to accommodate double precision variables. We have also made modifications to ensure that £, > 0 for all n by setting f,, = 10 s when f,, drops below 1() -s These parameterizations may need to change as more experience is gained with the procedure. We sometimes find that the P o l a k - R i b i e r e variant does not converge for the proposed new value for /.t. In this case we recycle our old value and attempt a new jump using Eq. (23). The program is set to carry out a specified number of walks on bt or to stop the walk when the least squares function L achieves a predetermined small value, typically on the order of

H.K. McDowell. A.M. Clogston / Chemical Physics 211 (1996) 91-100

10 -8 to 10-10. The program allows for interactively changing the values of 3 and the number of walks. At the point where /x has attained a value such that L is within our specified criterion, we fix tz and turn on /3 from 0 to 1 in small increments on the order of 0,01. We find that the final value of L is not strongly dependent on 13 when the Gaussian form in Eq. (21) is used. In general we anticipate that a more sophisticated walk on /, and /3 will be required. The time required to obtain a suitable imaging of p t ( w ) through f,, using an IBM 6000 model 320H workstation can be substantial and on the order of hours. At the present time we are searching for methods to enhance the convergence procedure for obtaining a suitable value of the least squares fimction L. In the following section we examine a model system to demonstrate the quality of the MaxEnt imaging procedure.

3. B e s s e l f u n c t i o n T C F

In previous work on linear vibrational chains [38], exact time correlation functions have been found of the form

0,(t) = 2J2( 2O2ot) /~o2t,

(24)

where J 2 ( x ) is a Bessel function and w o is a frequency. The exact spectral density p~(w) is given by

pl( w) =

"n'to~ 2~2 (

1-

=0,

95

2.5 -

1.0 loglo L

0.0-

~

0.5 -2.5IOgl0 L

I -5.00.0 -7.5.

-10.0-

.2.o

0'0

21o

41o

60

-0.5

Iogl0 P. Fig. 1. Behavior of least squares function in Eq. (14) and inlormation functional 1 in Eq. (13) with respect to regularization parameter p~ for Bessel function test ease.

constraint in Eq. (13). The flattening out of I to form a plateau region as /x increases is an indication that the MaxEnt process has driven p l ( w ) close to its exact form and therefore I nearly to its asymptotic value. In Fig. 2 we show the behavior of L and 1 for a separate stochastic walk to increase /x using the parameters above, but with /x driven to larger values. We find that numerical breakdown occurs. For this case breakdown is basically due to computational noise dominating the calculation of L to the point that L no longer changes and I decreases linearly as - ~ L .

t I/2 7;7 o ! , t°2

2.5

~o > 2o~ 0 .

i [3=O

(25)

1.0

O-

0.0

-2.5 -

-1.0

We thus obtain that cr~ " ' ) = [ 2 ( 2 m + l ) ! / m ! ( m + 2 ) ! ] W o " .

(26)

Using these moments, we carry out the MaxEnt imaging procedure described in Section 2. In Fig, 1 we show the behavior of L and I as a function of /x with M = 14, N = 80, co~ = 0.0, ws0 = 3.0(,o 0, a = 0.1, ftol = 10 8, and /3 = 0. Note from Eq. (13) that I = 1 f o r / , = 0. We see from the figure that the maximum value of I for fixed /x decreases as /, increases and the fit, as measured by L, improves. This is a direct result of imposing the

IOglO L

I -5 -

-7.5 -

-10

-2.0 [~=0

e

-3.0 I

-4.0

lOglO Fig. 2. Example of b r e a k d o w n of /J-function and information functional I as p a r a m e t e r / x increases in Bessel function test case.

96

H.K. McDowell, A.M. CIo,¢ston / Chemical Physics 2 l I (1996) 91 - I00

To understand the origin of the numerical instability in Fig. 2, we have carried out a stochastic walk to increase /x using the parameters above with the Simpson evaluation of integrals replaced with a 64point Gauss-Legendre quadrature. In Fig. 3 we show the behavior of L as a function of /z for this computational run in the region where the first signs of numerical breakdown occur. We compare these results with data from the computational run in Fig. 2. We find that at /x = 105 and L on the order of 10 -8 the two quadrature methods produce slightly different results. The difference is a measure of the numerical error in the computation of L. As we increase /x, we find that discontinuous jumps in L occur for both integration methods. The origin of these jumps is clear from Eq. (14). Basically, the quadrature methods are used to compute ~]/,,) and a difference taken with the exact value or(m). As /x increases, this difference eventually samples only the error in ~('~) due to the quadrature method. From Fig. 3 we deduce that this error occurs in about the fourth to fifth digit. Thus, as /x increases, it is quaranteed that L will at some value of /z become a " n o i s y " function sampling only the error in ~("). The jumps in Fig. 3 are simply a manifestation of this noise. Examining the computation of I in Eq. (18), we see that there is also an underlying error in I due to the chosen quadrature scheme. Furthermore, as L becomes noisy, it is mulitiplied by a large value o f / x in Eq. (18) which means that I becomes noisy. .8an

-8.1. " •

~n~

o

U

D Up

~t

U D~

D

[]

-8.2. Iogl0

ao

• ••~*D••~

o11. •

~

n

[]

U

0o•0

L -8.3-

-8,4-

-8.5

a •

Simpson

*~ • t

* *, •

Legendre

s.ls

sls

s.~s

log10 P. Fig. 3. Behavior of L-function at large /x for two integration procedures using Bessel function lest case.

1.5

--• ..........

!3=0 :."'"..

mOPl(m )

~\"""....... 0.5-

00

0

exact maxent guiding

_t___

/

015

i

1

i

1.5 ~/¢.00

i--

2

2.5

Fig. 4. Comparison of bath spectral density for Bessel function test case with M = 14 computed exactly, by MaxEnt imaging, and by the guiding distribution.

As a result, significant breakdowns as seen in Fig. 2 can occur. The numerical breakdown in our MaxEnt image procedure at large /x in a practical sense is not a significant problem since better quadrature methods can be employed when breakdown is detected. Furthermore, we find that the plateau behavior of I is found for all cases studied and serves as an indication of when a computation has converged. In Fig. 4 we plot pl(o)) as a function of co when M = 14, N = 8 0 , o) t = 0 . 0 , o)80=3.0o) 0, 6 = 0 . 1 , f t o l = 10 8, L = 9 . 1 2 0 × 10 -9 , / x = 3 . 1 3 7 × 105 , and /3 = 0. We point out that the MTGLE procedure currently used by us [7] agrees very closely with the exact curve in Fig. 4 and we do not plot this comparison. Note first of all the pronounced shift of the MaxEnt image from the guiding distribution. Note further that the MaxEnt image reproduces the sharp cutoff at w = 2 co0 of the exact spectral density in Eq. (25). The oscillatory behavior at the top of the peak is worth investigation given that the MaxEnt process should maximally remove information resulting from such behavior. To examine this behavior further, we turn on /3 to the value one to remove numerically the guiding distribution. The resulting spectral density is shown in Fig. 5. We find that the measure of the fit L improves to L = 2.303 × 10-9 but only small enhancements are made in the MaxEnt image relative to the exact result. To some extent, this is in keeping with the fact that the shape in Fig. 4 has already changed appreciably from the

97

H.K. McDowell, A.M. Clogston / Chemical Physics 211 (1996) 91-100

1.25

2.5

- 1.00, *

exact maxent~

M 14

0.02~ -2.5-

=

|ogl0 L

c°091(~) 0 . 7 5 - / ' ~ 0 . 5 0 0 . 2 5 _ .

-5.0-7.5-

0.00

0

0'5

~

~

15

j~

-10.0 -2.0

25

0'.0

~o/e)0 Fig. 5. Effect of removing guiding distribution from MaxEnt image using regularization parameter /3.

guiding distribution. Turning off the effect of the guiding distribution is not expected to be large. We have also carried out computations in which the two regularization parameters p, and /3 are turned on at the same time. We find essentially the same results as shown in Fig. 5. To further clarify this conclusion and to verify the MaxEnt image in Fig. 5, we have examined the dependence of the image on the parameters which drive the MaxEnt process. We first tested the tolerance parameter " f t o l " of the P o l a k - R i b i e r e variant. This parameter determines convergence of the variant procedure using the rule

211'-II

~< ftol( I 1' I + I 11 + ~ ) ,

(27)

where I' is a new computation of the information functional and e is a small positive number (10-L0 to 10-15). W h e n the condition in Eq. (27) is satisfied, the P o l a k - R i b i e r e variant is taken to be converged as described in Section 2. W e tested ftol in the range 10 -8 to 10 - t l and found that the results are converged on this parameter. We next tested the effect of increasing the number of Simpson rule points from N = 80 to N = 100 and found no effect. The regularization process is driven by the slow turn on of the parameter p.. We next examined the effect of changing the rate of change of /~ by reducing the jump parameter 8 in Eq. (23) from 0.1 to 0.05. Only small improvement in the spectral density p~(w) as regards smoothness were found. We turn finally to the effect of the number of bath moments M. In Fig. 6 we present the effect of reducing M from 14 to 9 on the behavior of L as a function of p- W e find that L is slightly reduced for

210 IOgl0 Ix

410

6.0

Fig. 6. Comparison for Bessel function test case of L-function with respect to parameter /x for two choices of number of moments M and /3 = 0.

a fixed value o f / x . On the other hand we find in Fig. 7 that I is not depressed nearly as much for M = 9 as M = 14. This can be viewed as a direct result of the number of constraints imposed on I through the requirement that M moments are fit to accuracy L. In Fig. 8 we examine the spectral density p~(w) for M = 9 at /3 = 1. The guiding distribution for M = 9 is included for comparison. W e see that the shape of pt(w) is changed from that in Fig. 5. Thus, we find that the principal parameter determining the shape of pt(w) is the number of bath moments M. We further conclude that the oscillations in p~(w) are intrinsic to the fit and are determined by the number of moments. It is also of interest to examine the behavior of the MaxEnt image of p l ( w ) as o J ~ 0 . From Eq. (17) we anticipate that p~(0) will not be zero when /3 = 1 since this requires /x to become large and positive; that is, the functional form of p l ( w ) as w - ~ 0 at

|.0M

14 0.5- - . . . . . . . . . . o ~ ~ = I

0.0-

-0.5

2.0

010

210 Iogl0 g

410

6.0

Fig. 7. Same as Fig. 6 for information functional I.

98

H.K. McDowell. A.M. Clogston/Chemical Physics 211 (1996) 91 100 1.25

In terms of the digitized, quadrature version we have

- -

exact

1.0o-I

N

maxent

O,(t) : ~, r,,of,, cos ~o,,f.

guiding

0.75.

6=1

m0Pl(m) 0.50. 0.25. / 0.00 m/m0

Fig. 8. Comparison of bath spectral density tbr Bessel function test case with M 9 computed exactly, by MaxEnt imaging, and by the guiding distribution.

/3 = 1 does not behave correctly. On the other hand we recognize that the imaging process is not intended to produce exact results but numerically satisfactory results as determined by the measure L. Thus, the digitized values of f , near w = 0 can reflect p,(0) being close to zero for /x sufficiently large in such a manner that p~(~o) is adequately represented. In our method we compute f , as we turn on /3. Thus, p,(O)=flg(O)=O from Eq. (21) and the results presented herein use this relation. On the other hand we note that the data points for p~(w) in Fig. 5 tend to lie above the true curve as w ~ 0. This behavior is apparently due to the inadequate ability of the functional form in Eq. (17) to exactly fit p l ( w ) near co= 0. By forcing p l ( 0 ) = 0 we can anticipate that the digitized curve for p](o)) will compensate in some manner such as displaying oscillatory behavior. To test this notion and to avoid a precise value of zero for pl(0), we have also carried out runs with g(0) = 10 - s , but find the same results. To some extent, this outcome is due to the fact that our imaging process, which is based on power moments, samples the high frequency end of p](co) and is less sensitive to the low frequency end. In future work we plan to include additional constraints based on reciprocal power moments such as in Eq. (22) to further examine the low frequency behavior. In M T G L E theory it is of interest to consider the behavior of 0t(t) which is given by ~c 01(f ) = f

do) a0

p,(oo) cos

cot.

(29)

n= 1

(28)

A plot of O~(t) using Eq. (29) and the spectral density i)l(w) in Fig. 5 with M = 14, N = 80, w~ 0.0, w s 0 = 3 . 0 w 0, (3=0.1, f t o l = 10 - s , L = 9 . 1 2 0 >( 10 -9, /Z = 3.137 × 105. and /3 = 1 is shown in Fig. 9 in comparison with the exact curve obtained from Eq. (24). The agreement is manifest and indicates that MaxEnt imaging can be used to image both the TCF and its spectral density. We conclude our study of the Bessel function TCF in Eq. (24) by considering the effect of using a different and better guiding distribution. W e know from our previous work that usually the time-packet M T G L E image is available for use in this regard. For the Bessel function case the M T G L E image is the same as the exact image and we, thus, feed in as the guiding distribution a digitized image of the exact spectral density in Eq. (25). The digitized image is an approximate image since we are only using 80 Simpson points. We find that the step of turning on /3 to unity causes a significant degradation in L. We compensate for this degradation by increasing /.t, after /3 is set to unity, until L achieves a desired value. The MaxEnt image for /3 = l with N = 80, M = 14, o) 1 = 0.0, ws0 = 3.0w0, (3 = 0.1, riot = 10 -8, L = 9 . 3 7 6 × l0 -7, and / x = 2 . 3 0 6 × 1 0 5 is shown in Fig. 10 and compared to a result using the Gaussian distribution g(w) in Eq. (21) with 6 = 0.05. The shape of the two curves is similar. We find that the Gaussian guiding distribution is more sensi-

1.5 1.0~ t/ 0.5

0

exact maxent

o

9= 1

1; COot

2;

30

Fig. 9. Comparison for Bessel function test case of exact TCF with MaxEnt image from Eq. (29) using pt(w) from Fig. 5.

H.K. McDowell, A.M. Clogston/Chemical Physics 211 (1996) 91 I00 1.25 - 1-

Gaussian guiding (.05)

*

~.

Exact guiding

= 0.75. ~O0Pl(~O) 0.5"

0.25.

03

015

iio

115

20

25

to/toO Fig. 10. C o m p a r i s o n of MaxEnt images c o m p u t e d using Gaussian (,5 = 0.05) and exact guiding distributions.

tive to the parameter 6 and requires a smaller value. Plots of 01(t) similar to the one shown in Fig. 9 for each of these cases show equally excellent fits to the exact O~(t) with differences appearing only at longer times where Of(t) has substantially decayed. We conclude from the overall shape of the spectral distributions in Fig. 10 that setting /3 = 1 removes the influence of the guiding distribution to within differences due to numerical errors arising from the parameters which control the turn on of /z and /3 and the specific path by which /x and /3 are turned on. Thus, we find that inclusion of the /3-term in the definition of 1 in Eq. (13) achieves the desired goal. On the other hand the results in Fig. 10 point to the need for determining an optimal process for turning on the parameters /x and /3. We have tested the MaxEnt image procedure described herein on the single-vibron system previously studied by us using the MTGLE imaging process [3] and find similar results to those reported herein for the Bessel TCF. Generally speaking, we find small oscillations in the fit of .p](w) to known results and excellent agreement for O](t).

4. S u m m a r y A key step in the quantum MTGLE approach is the determination of the bath TCF (~(t) and its spectral density P s ( w ) b y use of the expansion in Eq.

99

(10) in GLE functions. As shown in previous work [3,5,7], this expansion often produces regions in which the spectral density p~(w) is negative and oscillatory. This ill-tempered behavior can be controlled to some extent and in this regard we have introduced several methods to determine the control parameters o ) and o) in the GLE functions [3,5,7]. On the other hand the advent of a new procedure to invert the TCF moments to produce the bath spectral density is desirable and is satisfied by the MaxEnt imaging procedure described herein. In particular for a Bessel function test case [35] we have shown that MaxEnt imaging produces a satisfactory lineshape for pl(o)) with the image of 01(t) in excellent agreement with the exact TCF. Degradation occurs only at longer times where the amplitude of 0t(t) has significantly decreased toward zero. We have examined the effect of using the timepacket MTGLE image as an accurate guiding distribution. Within small numerical differences due to the effects of the various parameters values and the algorithms used to turn on the regularization parameters /x and ,(3, we find that effects due to the guiding distribution are removed for /3 = 1. "Ilaus, our strategy of using the information functional 1 in Eq. (13) to obtain an unbiased image of the spectral density p~(o)) is achieved while at the same time a guiding distribution can be used to aid in the nonlinear process of arriving at the best image. The quality of the images achieved herein indicate that the MaxEnt process described should be further explored and optimized. The MaxEnt method described herein is a new twist on the Gull-Daniell approach developed to allow the imaging of a spectral distribution from knowledge of its moments. The new features include fitting the known moments to within a tolerance criterion rather than the usual practice of fitting noisy data, and the use of a guiding function, which can include a prior distribution, to enhance numerical convergence of the nonlinear search process. We have furthermore introduced a numerical procedure for removing the guiding distribution from the final result. This removal process may have application to the standard Gull-Daniell approach. We conclude that MaxEnt imaging is an excellent procedure for constructing the MTGLE bath spectral density Of(w) from its known moments o-~~m).

1O0

H.K. McDowell, A.M. Clo~ston / Chemical Phy.~ics 2 l I (1996) 91 - 1O0

Acknowledgement We thank the Robert A. Welch Foundation, grants Y-1234 and Y-1303, for their support of this research.

References [1] B.J. Berne and G.D. Harp, Adv. Chem. Phys. 17 (1970) 63. [2] P.C. Martin, Measurements and correlation functions (Gordon and Breach, New York, 1968). [3] A.M. Clogston and H.K. McDowell, Phys. Rev. B 44 (1991) 4978. [4] H.K. McDowe[1 and A.M. Clogston, J. Stat. Phys. 67 (1992) 331. [5] H.K. McDowell and A.M. Clogston, Phys. Rev. B 46 (1992) 126. [6] H.K. McDowell, J. Chem. Phys. 86 (1987) 1497. [7] H.K. McDowell, P. Tsai and A.M. Clogston, Chem. Phys. Lett. 210 (1993) 274. [8] B. Buck and V.A. Macauley, eds., Maximum entropy in action (Oxford Science, Oxford, 1990). [9] R.G. Gordon, Adv. Mag. Reson. 3 (1968) 1. [10] R.G. Gordon, Adv. Chem. Phys. 15 (1969) 79. [11] S.F. Gull and G.J. Daniell, Nature 272 (1978) 686. [12] X. Zhuang, R.M. Haralick and Y. Zbao, IEEE Trans. Sig. Proc. 39 ( 1991 ) 1478. [13] A. Wragg and D.C. Dowson, IEEE Trans. Inform. Theory IT- 16 (1970) 226. [14] B.R. Frieden, J. Opt. Soc. Am. 62 (1972) 511.

[15] Y. Alhassid, N. Agmom and R.D. Levine, Chem. Phys. Lett. 53 (1978) 22. [16] M.R. Teague, J. Opt. Soc. Am. 70 (1980) 920. [17] R.D. Levine, J. Phys. A 13 (1980) 91. [18] R.D. Levine, Adv. Chem. Phys. 47 (1981) 239. [19] L.R. Mead and N. Papanicoleau, J. Math. Phys. 25 (1984) 2404. [20] N.S. Tzannes and P.A. Jonnard, Opt. Eng. 26 (1987) 1077. [21] L. Hagg and O. Goscinski, J. Phys. B 26 (1993) 2345. [22] A. Tagli,'mi, J. Math. Phys. 34 (1993) 326. [23] E.T. Jaynes, Proc. IEEE 70 (1982) 939. [24] B.R. Frieden, Proc. 1EEE 73 (1985) 1764. [25] J.E. Shore and R.W. Johnston, IEEE Trans. Inform. Thcory IT-26 (1980) 26. [26] R.W. Harrison, Acta Cryst. A 45 (1989) 4. [27] H.W. Engl and G. Landl, SIAM J. Numer. Anal. 30 (1993) 15O9. [28] R.N. Silver, D.S. Sivia and J.E. Gubernatis, Phys. Rev. B 41 (1990) 2380. [29] J.E. Gubematis, M. Jarrell, R.N. Silver and D.S. Sivia, Phys. Rev. B 44 (199I) 6011. [30] E. Gallicchio and B.J. Berne, J. Chem. Phys. 101 (1994) 9909. [31] N.J. Dusaussoy and I.E. Abdou, IEEE Trans. Sig. Proc. 39 (1991) 1164. [32] G. Bricogne, Acta Cryst. D 49 (1993) 37. [33] T. Oscarsson, J. Comput. Phys. 110 (1994) 221. [34] M. Abramowilz and 1.A. Stegun, eds., Handbook of mathematical functions (Dover, New York, 1968) p. 886. [35] S.R. White, Phys. Rev. B 44 (1991) 4670. [36] H.K. McDowell, J. Chem. Phys. 87 (1987) 4859. [37] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical recipes (Cambridge Univ. Press, London, 1986) section 10.6. [38] H.K. McDowell, J. Chem. Phys. 85 (1986) 6034.