sample-size hysteresis - American Chemical Society

6670

J . Phys. Chem. 1991, 95,6670-6675

have shown that, provided the energy of the flow is held fixed, the entropy lowering increases monotonically as the Reynolds number is increased, i.e. as the flow is driven further away from equilibrium. Thus, in this case the distance from equilibrium is correlated with the deviation of the entropy alone from its equilibrium value. When the Reynolds number exceeds a certain critical value, the flow changes from laminar to turbulent. It has also been found that, when a stable turbulent and an unstable laminar flow with the same Reynolds number above the critical value are compared, the entropy is larger and the rate of entropy roduction is smaller for the turbulent flow than for the laminar.' l9 It hence follows that both the entropy and the entropy production of the stable turbulent flow are closer to their equilibrium values than those of the corresponding unstable laminar flow. Summarizing the above discussion of the cubic autocatalator in an isothermal CSTR and of Poiseuille flows in a tube, one may

B1

(19) Klimontovich, Yu. L.; Engel-Herbert, H. Zh. Tekh. Fir. 1984, 54,

Downloaded by FLORIDA STATE UNIV on September 10, 2015 | http://pubs.acs.org Publication Date: August 1, 1991 | doi: 10.1021/j100170a054

440.

state the following systematic trends: (i) As the system is forced further away from equilibrium by changing a suitable control parameter, its entropy or entropy production or both always deviate more greatly from their values at equilibrium. (ii) If stable and unstable states coexist under given conditions, the stable state invariably has the entropy or entropy production or both that are nearer their equilibrium values than those in the unstable state. Such tendencies have also been found to prevail in several other types of open chemical systems, including nonisothermal ones; details will be reported elsewhere. Supplementary Material Available: Figure lS, the entropy difference LLS(the average entropy difference hs)between stationary (oscillatory) states and the corresponding equilibrium states ~ ; 2S,the entropy as a function of the residence time T ~ Figure production u in stationary states and the average entropy production a in oscillatory states as a function of the residence time 7rm (9 pages). Ordering information is given on any current masthead page.

Systematic Errors In Free Energy Perturbation Calculations Due to a Finlte Sample of Configuration Space: Sample-Slze Hysteresis Robert H. Wood,* Wolfgang C. F. Miihlbauer, Department of Chemistry and Biochemistry, University of Delaware, Newark, Delaware 1971 6

and Peter T. Thompson Department of Chemistry, Swarthmore College, Swarthmore, Pennsylvania 19081 (Received: December 3, 1990)

Although it is well known that the free energy perturbation procedure is exact when an infinite sample of configuration space is used, for finite sample size there is a systematic error resulting in hysteresis for forward and backward simulations. The qualitative behavior of this systematic error is first explored for a Gaussian distribution, then a first-order estimate of the error for any distribution is derived. To first order the error depends only on the fluctuations in the sample of potential energies, AE, and the sample size, n, but not on the magnitude of AE. The first-order estimate of the systematic sample-size error is used to compare the effciencies of various computing strategies. It is found that slow-growth, free energy perturbation calculations will always have lower errors from this source than window-growth, free energy perturbation calculations for the same computing effort. The systematic sample-size errors can be entirely eliminated by going to thermodynamic integration rather than free energy perturbation calculations. When AE is a very smooth function of the coupling parameter, A, thermodynamic integration with a relatively small number of windows is the recommended procedure because the time required for equilibration is reduced with a small number of windows. The present results give a method of estimating this sample-size hysteresis during the course of a slow-growth, free energy perturbation run. This is important because in these calculations time-lag and sample-size errors can cancel, so that separate methods of estimating and correcting for each are needed. When dynamically modified window procedures are used, it is recommended that the estimated sample-size error be kept constant, not that the magnitude of AE be kept constant. Tests on two systems showed a rather small sample-size hysteresis in slow-growth calculations except in the first stages of creating a particle, where both fluctuations and sample-size hysteresis are large.

I. Introduction There is a great need to estimate and control the errors introduced into free energy perturbation calculations by the various approximations necessary in implementing these calculations on modem computers.14 Random errors are easily estimated, while systematic errors are much more difficult to estimate since it is (1) Mezei, M.; Beveridge, D. L. Ann. N. Y.Acad. Sci. 1986, 182, 1 . (2) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (3) Pearlman, D. A,; Kollman, P. A. J . Chem. Phys. 1989, 90,2460. (4) Pearlman, D. A.; Kollman, P. A. J . Chem. Phys. 1989, 91, 7831.

necessary to understand the source of the error in order to estimate it. This paper is concerned with the systematic sample-size errors in free energy perturbation calculations produced when a finite sample of configuration space is used. In many calculations this systematic error leads to large hysteresis errors; that is, when the calculation is done in the forward and the reverse directions, the results are systematically different.*-' These systematic errors (5) Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Chem. Phys. 1989,129, 193. (6) Straatsma, R. P.;Berendsen, H. J. C. J . Chem. Phys. 1988,89,5876. (7) Wood, R. H. J . Phys. Chem., in press.

0022-3654191 12095-6670$02.50/0 0 1991 American Chemical Society

The Journal of Physical Chemistry, Vol. 95, No. 17, 1991 6671

Systematic Errors in Perturbation Calculations can be reduced by taking more time to do the calculations. However, if one understands the errors and can correct for their effects, the calculations can be made more efficient. In the present paper a distinction will be made between two kinds of hysteresis: time-lag and samplesize hysteresis. In order to define these two kinds of errors, it is necessary to describe the various procedures for calculating free energies. In window type free energy perturbation calculations, by either Monte Carlo (MC) or molecular dynamics (MD) simulations, the free energy of a given transformation is calculated as the sum of the free energies of various steps or window^".^-^ Thus

lo00

-

'

,001 1


i-0

where X is a coupling parameter that changes the potential from that of state a to state b as X goes from 0 to 1 in N steps or Ywindows"(A, = i/m, AE,* = (E,+l - E,)/RT, and AA is the Helmholtz free energy for the transformation from state a to b. The brackets ( )& indicate a canonical ensemble average generated by either the MC or MD techniques at Xi. In the "window growth" procedure after each change in A, neq steps of MC or MD equilibration are taken to allow the system to equilibrate to the new potential function, and then the ensemble average is computed by taking data for a given number of steps, ndl. Throughout this paper similar equations apply to AG calculated in an isothermal, isobaric ensemble. In the "slow-growth" window procedure the number of windows taken is very large, the number of equilibration steps can be set to zero, and the ensemble average for each window is approximated from a single configuration or time step (n,,, = 1). Under these circumstances eq 1 reduces to N- I

= is0 CAEi*

100

0

N- 1

AA* = M / ( R T ) = -Cln ( e ~ p [ - A E ~ * l ) ~ , (1)

AA*

I 10

Figure 1. The average hysteresis Hy for a sample of n independent AE's with a Gaussian distribution is plotted versus the standard deviation of the distribution, u. Symbols: 0 , n = 4; 0,n = 100 (the error bars are approximately the width of the symbols). The lines are the first-order estimate of the hysteresis, u2/2n for (from top to bottom) n = 1,4, and 100. Note that at low u the Gaussian distribution results approach the first-order estimates, while at high u the differences between sample sizes of 1, 4, and 100 are negligible.

systematic errors produced by finite samples of configuration space. 11. Systematic Sample-Size Errors for a Gaussian Distribution In this section the qualitative features of systematic samplesize errors are illustrated for a Gaussian distribution. This introduces the concepts needed and gives the correct qualitative behavior, even though the distributions normally found in molecular dynamics and Monte Carlo calculations of free energies are not Gaussian. To calculate the error for sample size n, we used a computer to select n independent AE* values randomly from a Gaussian distribution of AE* and calculate &4*, using

(2)

Because there are no equilibration steps, there is a systematic lag between the Hamiltonian and the system config~ration."~This lag gives rise to time-lag hysteresis for MD simulations or configuration-number-lag hysteresis for MC simulations. These were explored in a previous paper.7 It was shown there how (1) the lags could be estimated from the results of the simulation and (2) the free energy could be corrected for the errors introduced by these time lags.7 In the window growth procedure there is no time-lag hysteresis if the number of equilibration steps is large enough. In this case there are still well-known hysteresis effects that get larger as the window size (AX) in~reases.2~The origin of this hysteresis is that, as the window size increases, one needs larger sample sizes to adequately sample the distribution, because the exponential weighting that is used increases the importance of the low AE* side of the distribution.lg2 Figure 7.1 in Allen and Tildesley2 illustrates this effect. In the slow-growth procedure for free energy perturbation calculations, there is no exponential weighting (see eq 2), so one might expect that sample-size hysteresis is not present. In this paper we show that there is a systematic samplesize error for small samples that is present in both the slow-growth procedures and window calculations. The systematic error persists even when there is no apparent problem with sampling the AE*'s that contribute most to AA*. The error increases with decreasing sample size and is largest when a single sample is taken. Section I1 of this paper explores the qualitative features of this error when the AEi* values have a Gaussian distribution. The major result is that the errors are strong functions of the standard deviation but not of the mean of the distribution. Section 111 examines the cause of this systematic error and derives a first-order estimate of it for any distribution. Section IV gives some examples of the magnitude of the effect for typical systems. Section V shows that the same error is present when umbrella sampling is used. Section VI gives recommendations for controlling and reducing this error. There are, of course, additional errors due to finite sampling of configuration space that are random and thus easily estimated by standard methods.2 This paper is only concerned with the

n

&in* = -In [(l/n)Eexp(-AE,*)]

(3)

i= 1

This experiment is repeated m times, and the results are used to calculate the average value of AA* for the m trials, (AA.*>,. Equation 3 becomes exact as the sample size n goes to =,Z provided that there are no configurations for which E,* is infinite but E,+,* is finite.* To calculate the systematic error, we compare the average value of AA,* to the value for n = -, AA,*, which can be calculated exactly for a Gaussian distribution: AA,* = -In Jp"p(AE*) exp(-AE*) d(AE*) = p

- u2/2 (4)

where p(AE*) is a Gaussian distribution of AE's; i.e. p(AE*) = exp[-(AE*

-~ ) ~ / 2 u ~ ] / u f i

(5)

where p is the mean of AE* and u is the standard deviation. The systematic sample-size error is designated Hy because it produces sample-size hysteresis. Hy is estimated as the difference between the exact AA,* and the average value of AA,*. Hy a (AA,*), - AA,* (6) Table I and Figure 1 show the results of a variety of calculations of Hy for various values of n, m, u, and p. For these results the error depends only on u and not on p within the experimental error. This result is important because it allows dynamic control of the systematic sample-size error, and as we shall see below, it is true to first order for any distribution. The results in Table I and Figure 1 also show that Hy increases as the sample size decreases, and for small n it is large even when u is small. The error is largest for n = 1, since (M,),is then just p and Hy is u2/2 (see eq 4). The line for n = 1 in Figure 1 was calculated from this relation. The results show that for small n there are significant systematic errors, even when the AE's that contribute the most to AA seem to be adequately sampled. It should be remembered that n in eq 3 is the number of independent (8) Bennett, C . H.J. Comput. Phys. 1916, 22, 245.

Wood et al.

6672 The Journal of Physical Chemistry, Vol. 95, No. 17. 1991 TABLE I: Systematic Sample-Size Error,Hy, for a Sample of D Values of AE* Taka from a Gausshn Distribution n ma o Y Sawb HVC dI2d 100

4

3000 100

1000

1 2 2 2 4 9 9 30 90 1

2

2


1

1000

1000

6 9 90 0.5 2 6 6 2

0 0 30 -30 0 0 -30 0 0 0 0 0 0 0 0 0 0 30 0

0.63 0.37 0.37 0.37 0.16 0.06 0.06 0.032 0.023 0.79 0.71 0.58 0.55 0.50 1.00 1.00 1.00 1.00

0.013 (2) 0.20 (4) 0.13 (4) 0.15 (4) 2.0 ( I ) 22.6 (4) 21.9 (4) 382 (1) 3834 (4) 0.15 (2) 0.91 (4) 13.1 ( I ) 32.9 (2) 3956 (2) 0.06 ( I ) 1.4 (4) 15.4 (2) 15.2 (2) 1.97 (7)

0.005 0.02 0.02 0.02 0.08 0.405 0.405 4.5 40.5 0.125 0.5 4.5 10.1 1013 0.063 1 .00 9.00 9.00 2.00

error or hysteresis for a forward step. The equation for the error is m n

Hy = lim {In [(l/mn)CCexp(-AEg*)] -

In order simplify this expression we need the following result for n

-In [( 1/n)Cexp(-AEij*)] =

AE* values taken. from the distribution. Because of correlations many time steps or configurations are taken before an independent value of AE* is sampled in MD or MC simulations. Since it typically takes between 40 and 200 time steps to get an independent sample in an MD run,’ many window type free energy perturbation calculations have been made with n = 1 or 2, and no account has been taken of the systematic errors given by the present calculation.’ Table I also shows that Hy increases very strongly with u. A detailed examination of the results for large u showed that, as expected, the important parts of the distribution were not well sampled; i.e., only a few percent of the AE’s in the sample contributed Substantially to the value of PAn* calculated by eq 3. Under these circumstances a large amount of effort is put into producing the sample of AE’s, and only a few percent of them are really useful in the calculation. A rough estimate of this effect is given by the sampling efficiency, Self,defined as twice the number of AE’s in the sample that are lower than AAn* divided by the total number of AE‘s in the sample. For n > 10, Senis a rough but useful measure of the fraction of the sample that has contributed substantially to the value of PAn*. Table I gives the value of the average sampling efficiency for a variety of calculations. Clearly, runs with n > 10 in which the sampling efficiency is below 50% have relatively high Hy and, in addition, waste more than 50% of the computing effort. In view of the central limit theorem of statistics, one would expect that the distribution of AE‘s in a free energy calculation is more nearly Gaussian when there are large numbers of independent particles that contribute substantially to AE*. When the distribution is close to Gaussian, the methods of this section could be used to estimate the systematic sample-size error. How useful this estimate would be remains to be tested. 111. First-Order Calculation of Sample-Size Hysteresis In this section we use a similar procedure to calculate the systematic error in a sample of size n for any distribution. As

before, we draw n independent AE* values from the appropriate distribution and repeat this procedure m times. This gives values of AEij* for i = 1 to n and j = 1 to m. The systematic error in AA* is then calculated as the difference between the average of m experiments, each with a sample size n, and one experiment with a sample size mn. In the limit, as m m, the one experiment with sample size mn gives the exact result, while the m experiments with sample size n give the correct average result of using a sample size n. The difference between these two is just the systematic

-

(8s)

I= I

n

-In [ ( l / n ) Z ( l - AE,* i=1

+ AEij*2/2 - A E i j * 3 / 6 + ...)I

= (8b)

-In [ 1 - CAEij*/n

+ CAEij* ’/2n I

1

- CAEij* 3 / 6 n + ...I = 1

(8c)

x A E i j * / n - [CAEij*’/n- (CAEg*/n)2]/2 i

i

i

‘Number of independent sets of samples with each set containing n independent values. Average sampling efficiency as defined in the text for the m samples. CDefinedby eq 6. The number in parentheses is the standard deviation of the last digit of Hy estimated from the m values of Hy. “Predicted value of Hy from eq IOc. Note that this first-order prediction is accurate when either n or o is small.

j=li=l

P-

CAEij*/n - 8?(n - 1)/2n i

+ ... =

+ ...

(8d) @e)

where eq 8b is derived by expanding the exponential in eq 8a in a Taylor series. Eq 8d is derived by expanding the logarithm in eq 8c and collecting terms according to the power of Mu*.The second term in eq 8d is closely related to the usual unbiased estimate of the population variance,, !8 derived from the jth sample of n values of AE*. The definition of :8 is 5:

= C(AEi,* - (AE*ij))’/(n - 1) = i

[C(AEjj*)2 - n(AE*ij)’l/(n - 1) (9) i

To derive eq 8e we have substituted eq 9 into eq 8d. A similar expression holds for the single experiment with a sample size of mn and the unbiased estimate of population variance, a,,:. Substituting these two results into eq 7 gives

+ %,,,,,’(mn- 1)/2mn + C(l/m)[CAEij*/n - 8?(n - 1)/2n] + . .I j i

Hy = lim {-CCAEg*/mn nr-

i /

lim [amm,2(mn- 1)/2mn P-

- (8?)(n - 1)/2n + ...I u2/2n

+ ...

-

=

= (loa) (lob)

where (8); is the average value of :8 for the m samples. Equation 1Oc follows from lob, since in the limit m =, (8); and’ ,8 both go to the population variance, u2, even when the population is not Gaussian. Equation 1Oc gives the systematic sample-size error, Hy, for any distribution calculated to first order by Taylor series expansion. Thus, there is a systematic error that increases as n decreases and, to first order, depends only on the standard deviation of the sample and not on its mean. These are the same qualitative conclusions reached for a Gaussian distribution in section 11. We have used eq 10 to estimate the first-order error that we should observe for the experiments on a Gaussian distribution in Table I. The results (see Figure 1) show that Hy seems to approach the first-order estimate at low u as expected. For n I 4 and u I1, the first-order estimate (eq 1Oc) is accurate within experimental error, but as expected the fintsrder corrections are increasingly inaccurate as the standard deviation of the sample increases. The limiting behavior of In Hy as u becomes large is interesting. Figure 1 shows that for a Gaussian distribution when u is very large, the differences in In Hy for samples of different size (n) are small. As u increases In Hy appears to converge on the first-order value for n = 1, namely In (u2/2). For n = 1 there is a simple heuristic argument that illustrates the origin of the systematic sample-size error. In free energy


Systematic Errors in Perturbation Calculations perturbation calculations with large step size, the molecule or molecules being perturbed find themselves, on average, in an environment that is far from the equilibrium environment. This environment is almost always a higher energy environment than the equilibrium environment. Thus, most of the AE's calculated are higher than the AE a t equilibrium, so that the average AE is always more positive than the free energy. This is just the Gibbs-Bogoliubov Another way of looking at this is that the irreversible work of producing a perturbation is always larger than the reversible work, since the reversible work (and thus the correct free energy) is calculated in the limit of infinitely small step sizes. The irreverisibility of the process gets larger as the step size increases. The error of this irreversibility is reduced by exponential weighting, as in eq 1, and the use of larger sample sizes. However, finite samples do not eliminate all of the error. For the case where the standard deviation of the distribution is very large, it would be expected that the most important part of the distribution used in determining M,* is sampled very rarely. Under these circumstances the results of most calculations would yield a AA that is very high but occasionally an experiment would give a much lower value. Thus, the average error in AA would be quite different from either the median or the most probable error. In looking at the computer-generated distribution of results for a Gaussian distribution with a sample standard deviation of 6, it was found that the distribution of PA,* was skewed in the expected direction. We can use the present first-order estimate of the systematic sample-size error to calculate the error in averaging the results of a forward run and a backward run. For a given window in the forward direction, from Xi to A,+', the experimental value of is greater than the true value by $/2nt. Similarly, for the reverse window from Xi+' to Ai, the experimental result for L4,+l,i* = is on average greater than the true result by ~,~~/2n,+~.In the limit of small AX where the first-order correction is accurate, uiwill approach ui+',and the two hysteresis errors exactly cancel if the average of hA,,i+I for the two runs is taken and if n for the forward and reverse runs are equal for each window. This result shows that if n is kept the same for each window in the forward and reverse runs, then the average of the two runs is a much better estimate than either run by itself.

IV. Examples of Sample-Size Hysteresis Several short free energy perturbation calculations were run to assess the magnitude of the sample-size hysteresis. The first example was a slow-growth, free energy perturbation calculation of the mutation of half of the molecules in a box of 300 modified OPLS MeOH to TIP3P H20.I0 This perturbation was run with a modified AMBER program,"+12 which allowed multiple-molecule perturbations. The large number of molecules perturbed (150) leads to large AE, improved sampling, and reduced fluctuations. The perturbation was carried out from X = 0.99-0.96 in steps of AX = After each femtosecond time step in the molecular dynamics run, AE* was calculated for AX steps of f1W2, f l W 3 , *IO-', and *IF5.In this way the effect of step size on hysteresis a t exactly the same set of phase-space points was calculated. Time-lag hysteresis is not a factor because the time lag is the same for all of the calculations. The 3000 values of AE* for each of the step sizes were used to calculate the observed sample-size hysteresis by using Hy = ( P E S ) - PA* (11) where AA* is calculated by eq 3. We find in this calculation that (9) Smith, W. R. In Sturisricul Mechunics; Singer, K., Ed.;Royal Society of Chemistry: Letchworth, U.K., 1973. (10) MOhlbauer, W. C. F. Ph.D. Dissertation, University of Delaware, Newark, DE, 1990. (1 1) Singh, U. C.; Brown, F. K.; Bash, P. A.; Kollman. P. A. J. Am. Chem. Soc. 1987, 109, 1607. (12) Singh, U.C.; Weiner, P. K.; Caldwell, J. W.; Kollman, P. P. AMBER (UCSF), version 3.0. Department of Pharmaceutical Chemistry, University of California, San Francisco, 1986.

The Journal of Physical Chemistry, Vol. 95, No. 17, 1991 6673 the observed systematic sample-size error is Hy = 190(AE*)Ah. The estimated systematic sample-size error from eq IOc is Hy = 170(AE*)AX. Thus, the systematic sample-size error is small In either case, eq 1Oc is at AX = and large at AA = a good estimate of the observed systematic errors. It is well-known that particularly large fluctuations occur in the process of growing a nonpolar solute into water. The above theory predicts that with these large fluctuations the systematic errors will be large and, indeed, Pearlman and Kollman' found a very large hysteresis effect in the process of [ ] Xe(aq). Half of the total hysteresis occurred a t the very beginning of the calculation (0 < X < 0.2). The calculated AG's are larger than the average AG, which is the predicted trend for systematic sample-size error. It is important to note that this effect cannot be time-lag error because, for this run, the time-lag error would have the opposite effect.' Our short preliminary run for the same process showed that the hysteresis and the fluctuations were similar to the above MeOH H 2 0 example except a t low X. At X = At X = 0, Hy = 12000.01, Hy = 70(A,!?*)AX for AX I and Hy = 380(AE*)AX when 10" (AE*)AX when AX = I AX I IO-*. The systematic error is greater and no longer proportional to AX a t X = 0. Presumably, the high hysteresis is due to strongly repulsive configurations where the Lennard-Jones potential is highly nonlinear. Again, eq 1Oc gave good estimates of the systematic sample-size errors from the observed fluctuations. It seems likely that the well-known problem of making accurate calculations with slow-growth perturbation methods for this kind of process is due in part to sample-size hysteresis. Sample-size hysteresis, together with time-lag hysteresis, could be a possible explanation for a puzzling observation of Pearlman and Kollman, who found that in a slow-growth, free energy perturbation of [ ] CH,(aq), the calculated free energy change was 2.81 kcal mol-' with essentially no hysteresis, whereas other more accurate calculations gave AG as 2.0-2.4 kcal mol-' for the same modeL3 We postulate that in their slow-growth calculation, both systematic sample-size and time-lag errors were present and that they canceled, thus hiding the fact that the calculation was not very accurate. This possible cancellation of errors indicates the need for separate estimates of both time-lag' and systematic sample-size errors in slow-growth calculations.

-

-

-

V. Umbrella Weighting Torrie and V a l l e a ~ ' ~introduced J~ umbrella sampling to get around the problems that arise when the more negative values of AE* are not well sampled. Torrie and Valleau sampled from a weighted distribution with the weight, W,a function of A,!?*. They showed that under these circumstances AA* is given by

J

b

-In [ (exp(-AE* + Vw*)) w] + In (exp(Vw*)) (12) On the right side, we have defined Was exp(-Vw*) with Uw* = V w / k T . This shows that umbrella sampling is equivalent to exponential weighting of a different distribution. If the distributions of (-A,!?* Vw*)and Vw* have standard deviations different from zero, then the present calculations show that both terms will suffer from systematic errors when finite samples are used. In the general case it is not expected that the two contributions to systematic sample-size error would exactly cancel. Therefore, we conclude that there are systematic errors in umbrella sampling calculations when finite-sized samples are used and that these systematic errors could be estimated to the first order by examining the standard deviations of the distributions.

+

VI. Discussion and Recommendations Now that we have derived the first-order estimate of systematic sample-size error, we can explore various strategies for reducing this error and make recommendations as to the most efficient way (13) Torrie, G. M.; Valleau, J. P. Chem. Phys. Leu. 1974, 28, 578. (14) Torrie, G. M.; Valleau. J. P. J . Compur. Phys. 1977, 23, 187.

6614 The Journal of Physical Chemistry, Vol. 95, No. 17, 19‘91 to do free energy perturbation calculations. In the past, it has frequently been recommended that the step ~ size be reduced until AE* = AE/kT was less than 2 or 3 . ’ ~ ’ The present results show that to first order the systematic error does not depend on the mean value of AE* and is, rather, proportional to the square of the standard deviation of AE*. (A reviewer has pointed out that this can be easily verified by substituting AEi* = (AE*) 6Ei* in eq 3 and factoring out of the exponential to get AA,* = ( E * ) -In (exp(-sE,*)). This shows that (AE*)is an additive correction, and the remaining contributions represent only the variation in E, not its magnitude. Similarly, it can be shown that the same substitution in eq 7 shows that Hy is independent of (AEi*).) Thus, a better recommendation is to reduce the step size until the standard deviation u of AE is small enough to make the estimated systematic sample-size error (eq 1Oc) negligible. In Pearlman and Kollman’s dynamically modified windows procedure for controlling hysteresis errors, step size is continuously adjusted to keep the average AE* in each window cor~stant.~ The present results show that a better procedure is to keep u2/2nconstant since AE could increase with X while u2/2n decreases. This dependence on u2/2n explains the empirical observation that larger values of AE* are acceptable in some calculations, for example, in electrostatic charging calculations,16J7 whereas for the process of creating a hydrophobic particle in water only very small values of AE* are a~ceptable.~J’It is well-known that in the beginning stages of creating a particle, the fluctuations in the energies are very large, and this explains why extremely small step sizes are needed for this kind of cal~ulation.~ Since many existing programs are set up to do slow-growth and/or window free energy perturbation calculations, it is instructive to use the present results to compare the two calculations. For small values of AX, the first-order estimate of the systematic sample-size error will be accurate and we can compare the two procedures. Consider two calculations that require the same computer time. The first calculation (example a) is a free energy perturbation calculation with N windows with each window sampled for a time T. Let the standard deviation of the values of AE* in this calculation be independent of X and equal to ua, The second calculation takes about the same amount of computer time but uses 2N windows with sampling time T / 2 for each window and a standard deviation of AE* in this run ub = u a / 2 . The standard deviation is half as large because fluctuations are normally proportional to AX for small AE*. In the case where the sampling time Tis larger than the autocorrelation time, 7 , for an independent sample, we have the number of truly independent samples in run a, n, = T / T ,and the number in run b, nb = T / 2 r . We then calculate the total systematic sample-size error as the sum of the hysteresis errors for each window and find


+

Thus, the systematic error (in the limit that the first-order correction is correct) is independent of the number of windows. Of course, as the number of windows get small, the higher order corrections to the systematic sample-size error come into play and the error will get larger as N decreases. This is what is normally found in practice. In the case when Tis less than 7, we find that the larger number of windows produces less systematic sample-size error. In this case, although successive samples are highly correlated, the systematic sample-size error does not change because the effective value of n, = nb = 1. Under these circumstances, the systematic sample-size error decreases and is inversely proportional to the number of windows. As a result, the calculation with the least sample-size hysteresis is the slow-growth calculation where a window consists of only one M D step. Thus, if t is the time step used in the MD integration, a window growth procedure (IS) Postma, J. P. M.; Berendsen, H. J. C.; Haak, J. R. Furuduy Symp. Chem. SOC.1982. 17. 55. (16) Urban, J.’J.; Damewood, J. R., Jr. J . Chem. Soc., Chem. Commun., in press. (17) Jorgensen,W.L.;Blake, J. F.; Buckner, J. K. Chem. Phys. 1989,129, 193.

Wood et a]. with 7 > T has at least r / t more systematic samplesize error than a slow-growth procedure. Values of 7 / t of 50 to 100 are not uncommon.’ To complete this analysis, a consideration of both the random and time-lag errors is necessary. Consider first the random error. The total random error is calculated as the square root of the sum of the squares of the standard deviations of each window. The standard deviation of the average of each window is just the standard deviation of a single observation divided by the square root of n, the number of independent samples. For small values of the standard deviation it is found experimentally that the standard deviation of a single observation is normally proportional to AX. In this case the random error is independent of the number of windows. The time lag between the Hamiltonian and the system configuration also causes hysteresis’ because we are sampling a nonequilibrium configuration. These errors have been discussed in detail in a previous paper.’ In the present context the important question is whether one can reduce time-lag hysteresis by going to half as many windows and doubling the time allowed for equilibration, Tq. It is easy to see qualitatively the answer to this question in the two limiting cases. The first case is when Tcs is much less than the equilibration time, T ~so, that the percentage equilibration in each window is very small. In this case the time-lag errors accumulate over many windows and the accumulated time-lag error is essentially independent of the roughness of the change in X with time. As a specific example, compare a slow-growth simulation (nq = 0, nWll= 1) with one containing the same number of time steps, but with each window consisting of one equilibration and one data collection step (nq = 1, nmll= 1). This case is close to the limiting case where A changes continuously and the time-lag hysteresis error depends only on the total number of steps, and not the number of steps in each window. Figure 5a of ref 7 illustrates the case where the time-lag error accumulates over many time steps. In this case time-lag errors are independent of the number of equilibration steps. Therefore, it is better not to use any equilibration steps because the systematic sample-size error is smaller and the random errors are essentially unchanged because of the correlation between time steps. The other limiting case is when Togis much greater than rq. In this case it is easy to see that going to fewer windows greatly decreases the equilibration error. For example, if during Tw in a given window the equilibration error is reduced to 10% of its initial value, then going to half as many windows and, thus, doubling Tq will reduce the equilibration error to 1% of its initial value in each window. As far as equilibration errors are concerned this is a much more effective calculation and requires the same amount of computing time. Of course, for a small number of windows, AX is large and the systematic sample-size error is large so that free energy perturbation calculations done in this way are not recommended. However, as we shall see below, thermodynamic integration is the procedure of choice when a small number of windows is used. An alternate to the above procedures is to eliminate sample-size hysteresis entirely. The present calculation of the systematic sample-size error shows how this can be done. If one calculates values of AE/AX and takes the limit as AX 0, the systematic sample-size error vanishes. The result is just the well-known procedure of thermodynamic integration in which it is the derivative dE/dX that is calculated. There is no systematic sample-size error because the value of dE/dX sampled is the value at the X of the equilibrium ensemble. The disadvantage of this procedure is that the derivative can only be sampled a t a finite number of points, and so assumptions about the smoothness of the dE/dX as a function of X are necessary. It is not expected that the random errors in dE/dX will be very different from the random errors in AEIAX, so that random errors should not be affected by switching from free energy perturbation to thermodynamic integration. There is the advantage that systematic sample-size error is entirely eliminated, and this must be balanced against the disadvantage of requiring assumptions about the smoothness of dE/dX as a function of A. Of course, random errors

-

J . Phys. Chem. 1991, 95,6675-6683 and equilibration (time-lag) errors are still present. The finite difference thermodynamic integration of Meails avoids evaluating derivatives but introduces some sample size error. The result of this analysis, then, is that slow-growth free energy perturbation or thermodynamic integration, as well as thermodynamic integration with large step size, is recommended. It appears that slow-growth free energy perturbation and slow-growth thermodynamic integration will give the same results for many normal calculations, but for calculations where fluctuations are high and, in particular, for calculations involving the creation of a particle, the slow-growth thermodynamic integration is the safer procedure. In slow-growth free energy perturbation calculations, (18) Mezei, M. J . Chem. Phys. 1987,86, 7084.

6675

one should use fluctuations about a smooth curve through the results to estimate u and Hy. In this way one knows when systematic sample-size errors are affecting the results. Thermodynamic integration with large step sizes will be the best method when the functional dependence of AE on A is expected to be very smooth.

Acknowledgment. We thank Larry Hotchkiss for helpful discussions. The work was supported by the National Science Foundation under Grant No. CHE87-12204 and the Department of Energy under Grant No. DE-FG02-89ER14080. Such support does not constitute endorsement by N S F or DOE of the views expressed in this article. We thank Academic Computing and Instructional Technology at the University of Delaware-for support of the computations.


Molecular Dynamics Study of the Free Energy Functions for Electron-Transfer Reactlons at the Liquid-Liquid Interface Ilan Benjamin Department of Chemistry, University of California, Santa Cruz, California 95064 (Received: January 7 , 1991)

Free energy functions for electron-transfer reactions at the interface between immiscible polar and nonpolar solvents are investigated by using molecular dynamics simulation. The full free energy curves are calculated by using an umbrella sampling procedure for an electron transfer between two centers located either parallel or perpendicular to the interface, and these are compared with the free energy curves for the same process in the bulk polar solvent. In most cases,the parabolic reprentation of the free energy curve (obtained from its curvature near the minimum) overestimates the high free energy values, with the largest deviation observed for charged solute in the bulk. For the systems at the interface, the parabolic representation is somewhat better, and the curvature of the free energy curve near the minimum is shown to depend on the orientation with respect to the interface of the line joining the charge-transfer centers. The role of solute charge and solvation structure in determining the free energy curves for the different systems is discussed.

1. Introduction

The study of electron-transfer (ET) processes occurring at the interface between two immiscible liquids is of fundamental importance in chemistry, chemical engineering, and bi~logy.l-~For a full understanding of electron-transfer reactions at the interface, we need to consider many factors, some of which are applicable to general ET and some are unique to the interface. For example, for an ET to occur at the interface we must consider the diffusion of the donor and acceptor to the interface, their specific orientation due to interfacial forces, the change in the solvation structure between the bulk phase and the interface, and the change in the solvent polarization fluctuations. This is in addition to more general aspects such as the electronic structure of the reactants and their electronic coupling and its modification by the solvent and the role of intramolecular nuclear motion. The important role of the solvent in ET processes was recognized early in the development of the theories of ET reactions. In particular the concepts of solvent coordinate and solvent free energy functions play an important role in many experimental and theoretical investigations. These free energy functions were first calculated by Marcus (who built on the earlier work of Gurney' and Libbys) using a continuum dielectric approach6*'and ( I ) The Interfoce Structure and Electroehemicol Processes ut the Bounodry between Two Immiscible Liquids; Kazarinov, V. E., Ed.; Springer: Berlin, 1987. (2) Ulstrup, J. Charge Tronsfer Processes in Condensed Medio; Springer: Berlin. 1979. -. (3)'Gratzel, M. Heterogeneous Phoiochemicol Electron Transfer; CRC Press: Boca Raton, FL, 1989. (4) Gurney, R. W. Proc. R. Soc. London, Ser. A 1931, 134, 137. (5) Libby, W. F. J . Phys. Chem. 1952,56, 863. (6) Marcus, R. A. J . Chem. Phys. 1956, 24, 919. (7) Marcus, R . A. J . Chem. Phys. 1965, 43, 679.

0022-3654/91/2095-6675$02.50/0

statistical mechanics8 for ET in homogeneous solutions6 and for ET at metal-solution interface^,^.^ and recently for semiconductor-liquid and for liquid-liquid interfaces.I0J1 These expressions were used by Marcus to develop the theory of ET rate in solutions,12 at electrode^,^ and at liquid-liquid interfaces."J The earlier work was subsequently developed by Levich,I4 Dog~nadze,'~ and many other workers, using mainly the dielectric continuum formulation.* Although the macroscopic continuum approach has been very useful in explaining a large body of experimental data in a relatively simple way, it can be quite inaccurate, as has been demonstrated in many early studies of related systems and in recent experimental and theoretical work on solvation dynamic^.^^*'^ Therefore, a microscopic examination of the basic assumptions and the validity of the continuum approach is important. This can be done by using approximate integral equation theories,Is reduced dimensionality models,lg or the (potentially) exact computer simulation technique. The first application of computer simulation methods to understanding the free energy functions for electron transfer in bulk (8) Marcus, R. A. J . Chem. Phys. 1963, 39, 1734. (9) Marcus, R. A. J . Chem. Phys. 1963, 38, 1858. (10) Kharkats, Yu. I.; Vol'kov, A. G. J . Electroonol. Chem. 1985, 184, 435. (1 I ) Marcus, R. A. J . Phys. Chem. 1990, 94, 1050. (12) Marcus, R. A. J . Chem. Phys. 1956, 24, 966. (13) Marcus, R. A. J . Phys. Chem. 1990, 94, 4152. (14) Levich, V. G.Ado. Electrochem. Electrochem. Eng. 1966, 4, 249. ( I 5 ) Dogonadze, R. R. In Reactions of Molecules or Electrodes; Hush, N. S., Ed.; Wiley: London, 1971; p 135. (16) Simon,J. D.Arc. Chem. Res. 1988, 21, 128. (17) Bagchi, 8. Annu. Reu. Phys. Chem. 1989, 40. 115. (18) Calef, D. F.; Wolynes, P. G. J. Chem. Phys. 1983, 78, 470. (19) Hynes, J. T. J . Phys. Chem. 1986, 90, 3701.

0 1991 American Chemical Society

sample-size hysteresis - American Chemical Society

Recommend Documents