Article pubs.acs.org/jchemeduc
A Primer in Monte Carlo Integration Using Mathcad Chad E. Hoyer and Jeb S. Kegerreis* Department of Chemistry, Shippensburg University, Shippensburg, Pennsylvania 17257, United States S Supporting Information *
ABSTRACT: The essentials of Monte Carlo integration are presented for use in an upper-level physical chemistry setting. A Mathcad document that aids in the dissemination and utilization of this information is described and is available in the Supporting Information. A brief outline of Monte Carlo integration is given, along with ideas and pedagogy for optimum use of the Mathcad document in relation to application of this technique to quantum mechanical calculations to which undergraduate chemistry students are typically exposed.
KEYWORDS: Upper-Division Undergraduate, Physical Chemistry, Computer-Based Learning, Computational Chemistry
M
especially, students are usually familiar with the basics of numerical analysis (Riemann summation, Simpson’s rule, etc.). So, if students have seen these standard numerical integration techniques, then why is there a need for Monte Carlo integration? The answer to this is the well-known issue of the scaling of each technique. A simple exercise appropriate for an undergraduate student that explains how these standard integration techniques scale follows. To integrate via a Riemann summation, for example, one must take each coordinate (or axis) spanned by the function being integrated and discretize said coordinates into a number of points that your summation will “visit”. Consider the one-dimensional case where the x axis has been discretized into N points with the labels x1, ..., xN and the spacing between said points is dx. In this case, the integral I of the function f(x) is approximated by the formula
onte Carlo integration plays a key role in both equilibrium and dynamics calculations for chemically relevant systems.1−5 Thus, it becomes useful for undergraduate chemistry students to have at least a cursory understanding of the basics behind this technique. The general ideas behind Monte Carlo integration have widespread application and are not limited to integration,6−8 and a hands-on experience with this technique provides yet another useful and necessary mathematical skill to the student.9 This article provides a means to expose upper-level undergraduate students to a Monte Carlo framework using a Mathcad document, while concurrently providing a justification suitable for this audience. The context of the document (the calculation of quantum mechanical integrals that are essential to a physical chemistry curriculum) also allows for the reinforcement of concepts such as the Born interpretation of the wave function due to the probabilistic nature of the Monte Carlo technique. Integration, whether it is numerical or analytical, is obviously of primary importance when it comes to calculating probabilities relating to observables of quantum mechanical systems. From problems that predict the probability of finding a particle in a given position, to problems that ask for the average of some operator, expressions that include integration are numerous within the context of a physical chemistry curriculum. Providing skills that allow the students to analyze such expressions is, thus, essential. Although analytical calculations are commonplace, exposure to the procurement of numerical solutions is also important. Because of the inherent nature of quantum mechanics and the limitations that exist in regards to finding exact solutions, numerical solutions to complicated problems are oftentimes the only true means of approximating an answer. In the context of integration © XXXX American Chemical Society and Division of Chemical Education, Inc.
N
I=
∑ f (xi) dx i
(1)
As N is increased, the approximation becomes more and more valid, but each “visit” is more computational cost. If, for example, it took 10 points or “visits” to get an adequate answer for a one-dimensional system, what happens if one then decides to integrate a similar function that is two-dimensional? The extra dimension requires us to discretize 10 points on the x axis and the y axis both, which then produces a total of 100 points for the computer to “visit” to analyze the full behavior of the function. Such scaling is referred to as exponential scaling (102 = 100 by the above example; for 3 dimensions, 1000 points are needed as 103 = 1000, etc.), and this is the real drawback of this type of numerical integration. Every added dimension causes an
A
dx.doi.org/10.1021/ed400013d | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
distribution function in their analytical coursework, and even learn about Boltzmann distributions in the context of the kinetic molecular theory of gases. Monte Carlo allows us to build upon these already established skills, as the sampling function that we choose for a given integral must have the important properties indicative of probability distributions. Thus, it must be real and positive everywhere. The sampling function must also mimic the integral’s shape as closely as possible, as an ideal sampling function will be small where the function being integrated is small (thus making said points less important to the process) and the sampling function will be large where the function being integrated is large (thus making such points more important to the process). The sampling function, if being treated as a true probability distribution, must also contain 100% of the important parts of the function you are trying to integrate (and this is simply normalization). With this in mind, let us consider an example that will spell out the Monte Carlo process as a whole. Let I be an integral we are trying to evaluate, and let K(x) be the function we are integrating between limits xi and xf, then
exponential growth in computational cost. The severe impact of this problem can be illustrated if one considers the calculation of the average of an observable A for a system described by a ⎯r → ⎯ → ⎯ multidimensional wave function ψ(→ 1, r2 , ..., rn ) that depends on the coordinates of all n atoms of said system, ⟨A⟩ =
∫ d→⎯r1 ∫ d→⎯r2··· ∫ d→⎯rnψ *(→⎯r1, →⎯r2 , ···, →⎯rn)Â ⎯r , → ⎯ → ⎯ ψ (→ 1 r2 , ··· , rn)
(2)
If each ri⃗ is a position vector for the i-th particle of the system, then the particle position itself acts as a dimension in the integration. Even 1000 particles (not nearly a mol) would produce an integral that even well-equipped computers cannot handle in a timely fashion. This issue can also be addressed through the analysis of the error associated with these so-called “grid-based” methods that require discretization of the coordinates of integration. Such processes have an error that varies on the order of (1/N)r/d, where “d” is the dimension of the integral being evaluated, r is the number of bounded derivatives of the function being integrated, and N is the number of points the computer must “visit”.10 In this respect, again, we see that the true limitation prohibiting the application of such grid-based techniques to high-dimensional integrals lies in the fact that the error has an exponential dependence on the dimension. For example, let r = 1 and d = 2, then by this formula the error associated with N = 100 is 1/10 . If the dimension is increased by 2, to get the same error you would need to increase N by a factor of 100 as (1/ 10000)1/4 = 1/10. This severe impact on computational effort brought about, again, by an increase in dimension makes it extremely difficult if not impossible to apply standard numerical integration techniques to many integrals of interest in physical chemistry.
I=
∫x
xf
d x K (x )
(3)
i
Suppose we note a part of K(x) has the appropriate properties that would make it a good Monte Carlo sampling function. We then factor K(x) into a Monte Carlo sampling function, s(x), and everything leftover is the Monte Carlo integrand, m(x). I=
∫x
i
xf
d x K (x ) =
∫x
i
xf
d x s (x )m(x )
(4)
As an example of how one might do this, consider Figure 1, which illustrates the simple one-dimensional integral of the
■
METHOD Monte Carlo integration is used to overcome the deficiencies of conventional numerical analysis tools. Unlike the exponential scaling mentioned above, the error associated with performing a Monte Carlo integration is independent of the dimension, and in fact varies on the order of (1/N)1/2, thus making it more beneficial relative to the above-mentioned grid-based methods for integrals of high dimension.10 In this expression, N refers to the number of steps in the Monte Carlo process, which will be explained later, but note that the dimension of the integral being evaluated no longer plays a role in the error associated with the integration process. There are no “grids” in Monte Carlo integration and instead the process is built on the fact that any integral can be turned into one where certain statistical truths apply.11 The true advantage of Monte Carlo comes from realizing this and being able to take an integral and divide it up into two distinct parts: a sampling function and a Monte Carlo integrand. The process is carried out by choosing a sampling function that acts as a probability distribution that designates the importance of all the coordinates of the function being integrated. It is used to produce a set of said coordinates, weighted by how much they contribute to the overall integral, which are then fed into the Monte Carlo integrand for analysis. The average of the Monte Carlo integrand values produced from this set yields a numerical approximation for the integral of interest. Note that students are familiar with probability distributions, as they are exposed to the Born interpretation of the wave function, are introduced to the Gaussian probability
Figure 1. Example of the one-dimensional integral of the function x sin2(x) over a range of zero to π showing the breakdown of the functions within an integral based on choice of sampling function and Monte Carlo integrand.
function x sin2(x) over a range of zero to π. If one plots x sin2(x), one notices that the points near the end of the range contribute far less to the area under the curve when compared to the points near the center of the range. The points at the ends of the range are thus less important than the ones in the middle of the range to the overall integral, so we choose a sampling function that will enforce this truth. A plot of the components of the function being integrated shows that a choice of ‘x’ as our sampling function would not be ideal, as such a choice would make the points near π most important B
dx.doi.org/10.1021/ed400013d | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
even though the function being integrated is approaching zero here. However, one can see that sin2(x) mimics the shape of the total integrand closely (decreases at the ends of the range, is largest near the center), and thus becomes a natural choice for the sampling function. This choice leaves us with ‘x’ alone as our Monte Carlo integrand. Whatever we do in regards to manipulating this integral into a form appropriate for Monte Carlo integration, we have to be sure that we retain the equality to the original integral of interest. With this in mind, let us impose the fact that s(x) must be normalized. To do so, we introduce a normalization constant C to the integral I=
xf
∫x
I=
i
C C 1 = C
d x s (x )m (x )
∫x
xf
xf
dx[C·s(x)]m(x)
i
∫x
xf
d x C · s (x )
i
(5)
(6)
Having transformed our integral into one with the appropriate structure for Monte Carlo integration, we can go about using the sampling function to select a set of points that we will feed to our Monte Carlo integrand. Again, Monte Carlo is all about taking advantage of the statistics associated with a given integral. Thus, we treat s(x) as a probability distribution that determines the importance of the x coordinates over the range of interest, and so we want the important parts to be visited frequently. However, we also want the less important parts to be visited as well, albeit less frequently. To do so we use a method called importance sampling and introduce randomness into our analysis (again, taking advantage of the statistical nature of our imposed setup). If we are at a given x and we want to know the next x that will be fed into our Monte Carlo integrand, we take a step in a random direction and then analyze what is referred to as the “critical ratio”. Let us label our old coordinate as xold and our new coordinate as xnew, then the coordinates are related as xnew = xold + α(1 − 2·n1)
(7)
where n1 is a random number having a range of values between zero and one, and α is an adjustable step-size parameter chosen so that the % of steps rejected is near 50%, which is required for an accurate approximation. The critical ratio is then analyzed as the following inequality
s(xnew ) ≥ n2 s(xold)
1 N
N
∑ m(xi) i=1
(9)
IMPLEMENTATION With this simple outline, students with adequate proficiency in coding should be able to produce a computer program in the language of their choice that will integrate via Monte Carlo integration. To familiarize themselves with the process, simple one-dimensional integrals can be used to test such programs. Calculations involving quantum-mechanical averages become an ideal candidate for these tests due to the existence of |ψ|2 in the final integral expressions for these objects. This allows for a reinforcement of the Born interpretation of the wave function, as the |ψ|2 parts of these integrals are more than adequate choices for Monte Carlo sampling functions. Although being able to code the process is worthwhile, the true purpose of this article is to report on a Mathcad12−14 document, which is the result of an undergraduate research project, that performs Monte Carlo integration in the context of this mathematical software (available in the Supporting Information). Because of Mathcad’s prevalence15−17 and the straightforward programming tools embedded within its framework, it is a great platform for undergraduate chemistry majors to explore Monte Carlo integration. Students who are fluent with Mathcad could come up with such a document themselves if the instructor sees fit, but dissemination of the document presented here is beneficial due to its versatility and worth in regards to providing a simple introduction on this subject. It should also be noted that although this particular document is based on the use of Mathcad, other software packages are also adequate tools for the described exercise, such as Maple,18 Matlab,19 and Mathematica.20 To present and explain the Mathcad document, the essential parts of the code are broken into two distinct parts. Please note that all variables are defined with clarity in mind, and this clarity should afford easy translation to settings other than Mathcad. Box 1 shows the input part of the document, which defines the important parameters and functions necessary to complete the Monte Carlo process. The setup shown is for the example given in Figure 1. After setting the stepsize, upper and lower limits of integration and the variable NumOfTimes, which is the number of Monte Carlo steps that the document will take toward integration, one simply defines the integral of interest (for the purposes of comparison to an accepted value) along with a choice of sampling function and Monte Carlo integrand. Box 2 shows the main part of the Monte Carlo code. The Mathcad document uses the built in random number generator (rnd) to produce random numbers between zero and one where appropriate. The first position is chosen randomly and
and we determine C by the common integral 1=
d x K (x ) =
■
d x s (x )m (x )
i
∫x
xf
It is of interest to note that as N → ∞, the Monte Carlo solution will approach the exact answer. Although Monte Carlo does overcome some deficiencies related to Riemann summation, it would be negligent not to point out when Monte Carlo becomes unusable. For integrals that are very oscillatory, the use of Monte Carlo becomes problematic. This is the famous “sign problem” that prohibits the use of this technique in certain situations. However, many smoothing techniques and other approximations have been developed to help theoreticians continue to make the necessary calculations involved with understanding the quantum mechanics of chemical systems.
i
=
∫x
(8)
where n2 is another random number between zero and one. If this inequality is met, the move is accepted and xnew becomes a coordinate to be analyzed by our Monte Carlo integrand. If the inequality is not met, then the coordinate is not moved and xold is fed into our Monte Carlo integrand instead. To finalize our example, imagine we use our sampling function to select a set of coordinates in succession, each labeled xi where i = 1, 2, ..., N, N being the number of Monte Carlo evaluations or “points” in our process. The integral of interest is found by C
dx.doi.org/10.1021/ed400013d | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
array of information is produced. The first position in the array (designated as “Approximation” in the code) gives the Monte Carlo approximation to the integral of interest. The second portion of the output array (designated “RejCounter”) is the number of rejected steps that can be analyzed to ensure that the process produces a rejection ratio between 40% and 60%. The third component of the array (designated ‘Xn’) contains the series of points chosen by the Monte Carlo sampling function (oftentimes referred to as the “random walk”) that can be plotted for a visual insight in regards to how the sampling function chooses the points that will be fed into the Monte Carlo integrand. The output is shown as {991, 1} due to the fact that there were 990 total steps in the random walk for the example shown, and each step is accounted for in the output. The fourth position in the array (designated “Approx”) is the step-by-step account of the accumulation of the final integral result via the 990 Monte Carlo passes used in the given example. This last component can be coupled with the random walk to provide a detailed picture of just how the final integration result is produced.
Box 1. The input section of the Mathcad document performing the Monte Carlo integration
Box 2. The main part of the Mathcad document, which performs the Monte Carlo integration. Output is in the form of an array of information shown in the top right corner.
■
STUDENT EXERCISES To verify the validity of this Mathcad document, a group of third-year chemistry majors in their second semester of physical chemistry were asked to perform a series of tasks relating to its use. During a three-hour session, students were introduced to the basic ideas of Monte Carlo, and then the Mathcad document shown in Box 1 and Box 2 was presented to them. As the code takes only a few seconds to run, the students were first instructed to analyze a printed version of the code and annotate it to include the purpose of each line and how it relates to the Monte Carlo introduction they were exposed to. Once this was accomplished, a set of integrals was presented to the students to complete, such as those given in Table 1. A simple quadratic integral was done first, followed by integrals that the students were initially exposed to when presented with the particle in a box and harmonic oscillator quantum mechanical systems. For each integral, the students started by choosing a sampling function and Monte Carlo integrand. They then proceeded by setting a stepsize and analyzing the percent of rejected steps. If the percent rejected steps was too small, they increased the stepsize, and if it was too big, they lowered the stepsize. By doing this, in a trial and error fashion, the students were able to achieve the required rejection rate. They then performed convergence tests by varying the number of Monte Carlo passes. They quantified their error in the form of a percent error analysis for single runs, and then they also tested the stability of their results by seeing how the error changed with the use of different random number sequences. The multiple results for a single integral were then averaged and a standard deviation was found (which is the common practice for integrals where the correct solution to the integration is unknown). Once this was complete, the students were asked to prepare a formal lab report outside of class that included their
the main Monte Carlo sequence is initiated. Steps are taken and the new position is first tested to see if it falls within the range of the integral. If it does not, the move is automatically rejected through assigning the position a value at the minimum of the sampling function (in this case, zero). The critical ratio is evaluated, and the approximation for the integral is calculated in the fashion explained above. Once the sequence is complete, an
Table 1. Student Results of Using the Mathcad Document for a Selection of Integrals Integral
s(x)
m(x)
Num of Times
Stepsize
Rejected Steps (%)
Monte Carlo Result
Percent Error
∫ 40 x2 dx ∫ 50 x[sin(πx/5)]2 dx −(x−2)2 ∫ 10 dx −10 x e 2 2 −(x−2) 10 ∫ −10 x e dx
x [sin(πx/5)]2 2 e−(x−2) 2 −(x−2) e
x x x x2
990 990 990 990
4.1 5.2 4.5 5.5
44.24 45.25 52.32 57.47
21.567 6.229 3.584 7.982
1.100 0.334 1.102 0.075
D
dx.doi.org/10.1021/ed400013d | J. Chem. Educ. XXXX, XXX, XXX−XXX
Journal of Chemical Education
Article
(4) Nakayama, A.; Makri, N. Symmeterized Correlation Functions for Liquid para-Hydrogen using Complex-Time Pair-product Propagators. J. Chem. Phys. 2006, 125, 024503. (5) Kegerreis, J.; Nakayama, A.; Makri, N. Complex-Time Velocity Autocorrelation Functions for Lennard-Jones Fluids with Quantum Pair-Product Propagators. J. Chem. Phys. 2008, 128, 184509. (6) Bentenitis, N. A Convenient Tool for the Stochastic Simulation of Reaction Mechanisms. J. Chem. Educ. 2008, 85 (8), 1146. (7) Rabinovitch, B. The Monte Carlo method: Plotting the Course of Complex Reactions. J. Chem. Educ. 1969, 46 (5), 262. (8) Davis, S. A Variational Monte Carlo Approach to Atomic Structure. J. Chem. Educ. 2007, 84 (4), 711. (9) Gasyna, Z.; Rice, S. Computational Chemistry in the Undergraduate Chemistry Curriculum: Development of a Comprehensive Course Formula. J. Chem. Educ. 1999, 76 (7), 1023. (10) Gerstner, T.; Griebel, M. Dimension-Adaptive Tensor-Product Quadrature. Computing 2003, 71 (1), 65−87. (11) Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087. (12) Cady, M. P.; Trapp, C. A. A Mathcad Primer for Physical Chemistry; W. H. Freeman and Co.: New York, 2000. (13) Noggle, J. H. Physical Chemistry Using Mathcad; Pike Creek Publishing: Newark, DE, 1997. (14) PTC Home Page. www.ptc.com/product/mathcad (accessed Jun 2013). (15) Zielinski, T. Helping Students Learn Mathematically Intensive Aspects of Chemistry. J. Chem. Educ. 2004, 81 (1), 155. (16) Zielinski, T. SymMath, Fair Use, and Additions to the Collection. J. Chem. Educ. 2005, 82 (1), 172. (17) Krenos, J. Physical Chemistry: Thermodynamics (Horia Metiu); Physical Chemistry: Statistical Mechanics (Horia Metiu); Physical Chemistry: Kinetics (Horia Metiu); Physical Chemistry: Quantum Mechanics (Horia Metiu). J. Chem. Educ. 2008, 85 (2), 206. (18) Maplesoft Home Page. www.maplesoft.com/products/maple (accessed Jun 2013). (19) Mathworks Home Page. www.mathworks.com/products/matlab (accessed Jun 2013). (20) Wolfram Home Page. www.wolfram.com/mathematica (accessed Jun 2013).
own explanation of how Monte Carlo integration works while also reporting their integration results with an argument regarding why their personal choice of the sampling function for each integral was appropriate. Justifying their choice of sampling function leads naturally to a discussion of the Born interpretation as the |ψ|2 terms present in each quantum mechanical integral ultimately outperform any other possible choices. A set of typical student results is shown in Table 1. It should be pointed out that the number of Monte Carlo passes has a maximum value of 990 for this document due to an inherent inability for Mathcad to handle large numbers of iterations through the established loop, but even with this relatively low number of points, the percent errors generated are quite low. The results shown are single runs and the percent error is based off of comparing the Monte Carlo result to the “Accepted” variable (see Box 1) calculated using Mathcad’s integration tool.
■
SUMMARY This article and the Mathcad document (available in the Supporting Information) act as a primer and the outline of the document use that was presented in the article provides a foundation for further investigation. To fully appreciate the worth of Monte Carlo integration, those students with advanced coding skills can use this as a foundation to tackle multidimensional integration that follows the exact same rules as those presented (once a multidimensional integral is chosen and factored into a sampling function and Monte Carlo integrand in Box 1, minor changes are made in Box 2, which include taking random steps in each dimension of integration while enforcing the boundary conditions for the new dimensions set therein as well). Comparison of the computational cost of doing Riemann summation integration relative to Monte Carlo integration can be gleaned and the true worth of the process can be established. However, as a means to introduce and establish this interesting and widely used mathematical tool, the document is more than sufficient.
■
ASSOCIATED CONTENT
S Supporting Information *
Mathcad document. This material is available via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Ceperley, D. M.; Mitas, L. Quantum Monte Carlo Methods in Chemistry. In New Methods in Computational Quantum Mechanics; Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; Wiley: New York, 1996; Vol. XCIII. (2) Makri, N. Forward-Backward Quantum Dynamics for Time Correlation Functions. J. Phys. Chem. 2004, 108, 806−812. (3) Makri, N.; Nakayama, A.; Wright, N. Forward-Backward Semiclassical Simulations of Dynamical Processes in Liquids. J. Theor. Comp. Chem. 2004, 3, 391−417. E
dx.doi.org/10.1021/ed400013d | J. Chem. Educ. XXXX, XXX, XXX−XXX
