Journal of Chemical Physics
Feedback | Help | Exit
The Journal of Chemical Physics
   
 
 
 
Previous Article
On the rate of absorption of Brownian particles by a black sphere: The connection between the Fokker–Planck equation and the diffusion equation
It is shown that, if one takes the first two moments of the Fokker–Planck equation, one obtains, in the lowest order, the telegraph equation irrespective of whether one uses a full-range or a ha...
Next Article
The effect of frequency dependent friction on isomerization dynamics in solution
We investigate the effect of the frequency dependence of friction on the rate of photochemical isomerization in solution. Recent experiments have shown a significantly smaller effect of viscosity on r...

Extended Brownian dynamics. II. Reactive, nonlinear diffusion

J. Chem. Phys. 78, 2713 (1983); doi:10.1063/1.445002

Issue Date: 1 March 1983

You are not logged in to this journal. Log in

Gene Lamm and Klaus Schulten
Physik Department, Technische Universität München, 8046 Garching and Max Planck Institut für Biophysikalische Chemie, 3400 Göttingen, Federal Republic of Germany
An improved version of the ``extended Brownian dynamics'' algorithm recently proposed by the authors [J. Chem. Phys. 75, 365 (1981)] is given. This Monte Carlo procedure for solving the one-dimensional Smoluchowski diffusion equation is statistically exact near a boundary for a constant force and approximately correct for a linear force. The improved algorithm is both more accurate and simpler than the earlier version. In addition, the algorithm is extended to include diffusion near a reactive boundary or in a reactive optical potential. The treatment of diffusion for nonlinear forces is conveniently handled by choosing the time for a single diffusive jump locally. The algorithm converges as this jump time approaches zero. The appropriate modifications necessary to treat diffusion between two (possibly reactive) boundaries or diffusion with a spatially varying diffusion coefficient are also given. Finally, it is shown how multidimensional diffusion in a spherically symmetric force field may be treated by the one-dimensional algorithm described here. As in the earlier paper, numerical results are presented and compared with analytical and numerical descriptions of the diffusion process to demonstrate the validity of the algorithm. The Journal of Chemical Physics is copyrighted by The American Institute of Physics.
History: Received 10 September 1982; accepted 26 October 1982
Permalink: http://link.aip.org/link/?JCPSA6/78/2713/1
BUY THIS ARTICLE   (US$28)
Download PDF (1915 kB) View Cart

KEYWORDS and PACS

Keywords
PACS
  • 05.40.+j
    Statistical physics and thermodynamics Fluctuation phenomena, random processes, and Brownian motion
  • 05.60.+w
    Statistical physics and thermodynamics Transport processes: theory
  • YEAR: 1983

PUBLICATION DATA

ISSN:
0021-9606 (print)   1089-7690 (online)
Publisher:
AIP is a member of CrossRef AIP

REFERENCES (40)
View in Separate Window/Tab
For access to fully linked references, you need to log in. For access to fully linked references, you need to Log in.
  1. D. L. Ermak, J. Chem. Phys. 62, 4189, 4197 (1975);
    2. D. L. Ermak and J. A. McCammon, ibid. 69, 1352 (1978).
  2. K. Schulten and I. R. Epstein, J. Chem. Phys. 71, 309 (1979).
  3. G. Lamm and K. Schulten, J. Chem. Phys. 75, 365 (1981); this paper is referred to as I throughout the present article and errata are listed in Appendix E of the present paper.
  4. J. Crank and P. Nicolson, Proc. Cambridge Philos. Soc. 43, 50 (1947);
    5. J. Crank, The Mathematics of Diffusion, 2nd ed. (Clarendon, Oxford, 1975), p. 141; see also Ref. 5 below.
  5. Z. Schulten and K. Schulten, J. Chem. Phys. 66, 4616 (1977).
  6. Throughout this paper the force will be expressed in units of kBT.
  7. M. V. Smoluchowski, Phys. Z. 17, 557, 585 (1916); see Eqs. (52)–(54).
  8. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (Dover, New York, 1972), p. 297.
  9. If b<0, distribution p2 becomes negative and two modifications of the procedure are necessary. First, in the partitioning phase, the random number r must be chosen from the interval (0,N), where N = N0+N1+|N2|. Second, the contribution due to endpoints distributed according to p2 (which will be small) must be subtracted from the total distribution of end-points. These complications are minor and are easily accounted for. This case occurs when dealing with a nonlinear force (Sec. VII) or with two boundaries (Sec. VIII). Appendix B gives a more general analysis of the problem of nonrectangular sampling.
  10. J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods (Wiley, New York, 1964), p. 36. This is the standard reference in Monte Carlo work; for a more recent summary see Ref. 35 below.
  11. We have programmed this algorithm and the results agree in all details with those presented here. Although It is slightly easier to program and more accurate (though not noticeably so), the algorithm is about 30% slower than the procedure described below.
  12. We have not observed any error near the boundary, even for extremely large b, as was the case with the algorithm in paper I. It is also possible to correct for the approximation by assigning a weight p2(x,t|x0)/[a exp(−b*x)] to the jump (see Appendix B), but the error is so small that this is unnecessary.
  13. Note the corrections to the figure of paper I given in Appendix E of the present paper.
  14. An application of the algorithm in which it is advantageous not to bias motion in the direction of the force is given in Sec. VII.
  15. We have also programmed this algorithm and have obtained results identical to those shown here. We regard this procedure as unsatisfactory, since for k>b/2 distribution Eq. (3.5) becomes negative and the modifications given in Ref. 9 above become nesessary. When k is relatively large, p2 must be accurately described in order to avoid negative distribution values for large x owing to the stochastic nature of the algorithm. In such cases, we have found Eq. (2.15) inadequate and have had to use Eq. (2.18) to represent the error function. This, in turn, requires additional partitioning with the result that the algorithm is about 30% slower and more difficult to program than that described here. An alternative procedure is to include the weighting factor mentioned in Ref. 12 but this then suggests what follows.
  16. In Ref. 5 the difference approach is applied to reactive diffusion.
  17. G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 (1930).
    18. [This paper is reprinted in Selected Papers on Noise and Stochastic Processes, edited by N. Wax (Dover, New York, 1954).] See also Ref. 7, Eq. (51).
  18. A rough, semiempirical restriction on the jump time is ct<ln(1+cx0/b), for b,c>0; see also Sec. V.
  19. In this section only, p0(x,t|x0) denotes some exactly known (unperturbed) distribution, and not the boundary-free distribution of Eq. (2.5).
  20. A. Szabo, K. Schulten, and Z. Schulten, J. Chem. Phys. 72, 4350 (1980);
    21. see also P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill, New York, 1953), Vol. I, p. 870.
  21. The diffusion domain is assumed to be (0,[infinity]). The surface terms at x = 0 vanish if p and p0 obey identical boundary conditions.
  22. L. I. Schiff, Quantum Mechanics, 3rd ed. (McGraw-Hill, New York, 1968), p. 304, Eq. (36.18).
  23. Note that, although the spatial factor zeta(x,x0) depends on the potential difference between the initial point x0 and the trajectory endpoint x, it must be computed at all intermediate jump points since the local force F0(x) is position dependent.
  24. Reference 10, pp. 57–59.
  25. A good introduction to path integrals is D. Falkoff, Ann. Phys. N. Y. 4, 325 (1958).
  26. W. Feller, An Introduction to Probability Theory and Its Applications, 3rd ed. (Wiley, New York, 1970), Vol. I, p. 470, Vol. II, p. 322, Eq. (1.3).
  27. In fact previous MC algorithms were unable to treat a reactive boundary correctly and an optical potential was required to describe reactive diffusion (Ref. 2).
  28. N. S. Goel and N. Richter-Dyn, Stochastic Models in Biology (Academic, New York, 1974), pp. 251.
  29. Equation (8.14) can also be derived assuming multiple reflections of the diffusing particle off the boundaries. This explains why Eq. (8.13) “just happens to be” normalized with these lambda values.
  30. When a force is present b at one of the boundaries will probably be negative. This requires the slight modifications to the algorithm given previously in Ref. 9.
  31. Although this algorithm has not been programmed, we estimate the additional complications would slow it down by 10%–30%.
  32. Reference 10, pp. 29–31. A general discussion of random number generation may be found in T. E. Hull and A. R. Dobell, SIAM Rev. 4, 230 (1962)
    33. and in B. Jansson, Random Number Generators (Victor Pettersons Bokindustri Aktiebolag, Stockholm, 1966).
  33. P. A. W. Lewis, A. S. Goodman, and J. M. Miller, IBM Syst. J. 2, 136 (1969). Additional pseudorandom number generators may be found in Ref. 8, p. 950.
  34. This particular algorithm may be machine dependent.
  35. F. James, Rep. Prog. Phys. 43, 1145 (1980) gives a recent and particularly relevant discussion; see pp. 1173–1175.
  36. C. Hastings, Approximations for Digital Computers (Princeton University, Princeton, 1955), pp. 169; this approximation is also quoted in Ref. 8, p. 299, formula (7.1.26).
  37. Reference 36, p. 192.
  38. M. E. Muller, J. Assoc. Comp. Mach. 6, 376 (1959).
  39. Reference 26, Vol. I, p. 35, 148.
  40. In Ref. 2, the rates given in Table I were calculated in this way for both the analytical and the Monte Carlo entries, i.e., the analytical entries were determined by employing in Eq. (D.3) the analytical yields from Eq. (4.13) of Ref. 5 rather than from the analytical expression (4.14) of Ref. 5 for −dN/dt as stated.

CITING ARTICLES
For access to citing articles, you need to log in.
For access to citing articles, you need to Log in.


Copyright © American Institute of Physics