Using the Metropolis Algorithm To Calculate Thermodynamic Quantities

Mar 8, 2011 - the Metropolis algorithm can be used. The method for doing this forms the main part of this article. If a thermodynamic “system” is ...
0 downloads 0 Views 1MB Size
Download PDF
ARTICLE pubs.acs.org/jchemeduc

Using the Metropolis Algorithm To Calculate Thermodynamic Quantities: An Undergraduate Computational Experiment Godfrey S. Beddard* School of Chemistry, University of Leeds, Leeds LS2 9JT, United Kingdom

bS Supporting Information ABSTRACT: Thermodynamic quantities such as the average energy, heat capacity, and entropy are calculated using a Monte Carlo method based on the Metropolis algorithm. This method is illustrated with reference to the harmonic oscillator but is particularly useful when the partition function cannot be evaluated; an example using a one-dimensional spin system is examined. The calculations are performed entirely using a computer algebra package, and the calculations are straightforward enough for undergraduates to investigate for themselves in the laboratory. KEYWORDS: Upper-Division Undergraduate, Physical Chemistry, Computer-Based Learning, Statistical Mechanics, Theoretical Chemistry, Molecular Modeling N, V, T) is defined as the “sum over states”

I

n classical thermodynamics, heat capacity is defined either as CV = (∂U/∂T)V or as Cp = (∂H/∂T)p. Using a calorimeter to measure the heat capacity of a substance is a common experiment in many undergraduate teaching labs. In other teaching laboratory experiments, students measure various spectroscopic constants and calculate partition functions, thus obtaining the vibrational or rotational contribution to a molecule’s heat capacity. The internal energy U at a given temperature and constant volume can be calculated from the canonical partition function Z as   2 D lnðZÞ ð1Þ U ¼ U0 þ kB T DT V

Z¼

U - U0 þ kB lnðZÞ T

ð2Þ

There is, however, another way of calculating internal energy and heat capacity, and this is via the relationships ÆEæ  U - U0, and CV ¼

1 ðÆE2 æ - ÆEæ2 Þ kB T 2

ð3Þ

where ÆEæ is the average energy and ÆE2æ is the average of the square of the energy.2 If the partition function can be calculated, then the averages in eq 3 are not needed. In many cases, however, there are obstacles to calculating the partition function and then the Metropolis algorithm can be used. The method for doing this forms the main part of this article. If a thermodynamic “system” is in contact with a heat bath, then the canonical (ensemble) partition function (at constant Copyright r 2011 American Chemical Society and Division of Chemical Education, Inc.

n

B

ð4Þ

where the energy of level n is En, its degeneracy is gn and kB is Boltzmann’s constant. For simplicity it is now assumed that gn = 1. The summation forming the partition function covers all n values, and even if n f ¥, the sum can often be evaluated algebraically. However, even if no algebraic solution is forthcoming, as in the case of the anharmonic (Morse) oscillator, Z can be evaluated numerically because with increasing energy the exponential terms very rapidly become insignificant and hardly contribute to the sum. In practice, therefore, the partition function summation has a limited rather than infinite number of terms. To appreciate this, if a harmonic oscillator has a frequency hν = 100 cm-1 (v = cν h), comparing terms with n = 0 versus n = 15 produces e-hν/(2kBT)/e-31hν/(2kBT) ≈ 1330 at T = 300 K. This ratio increases dramatically as ν increases or T decreases. Next, consider the more involved case of calculating the partition function for the hindered rotation about bonds in an alkane that might be part of a polymer chain or between adjacent residues in a protein. Imagine a polymer with 100 repeating units. Suppose that steric interaction between groups of atoms on adjacent units produces three potential energy wells and that there is only one energy level in each well. At any given time, the polymer chain can potentially exist in any one of 3100 ≈ 1047 configurations. The partition function cannot be estimated by summation until all terms are added, because the sum does not tend to a limit after only a few terms but is complete only when all are included. To put 1047 configurations into context, there have only been ≈4  1017 s since the universe was created, and even if it were possible with a supercomputer to calculate 1015 configurations each second, this single calculation

where U0 is the internal energy at T = 0. The heatRcapacity CV can be calculated by differentiating U; then using S = (CV/T) dT the entropy is obtained1 viz., S¼

∑n gne-E =ðk TÞ

Published: March 08, 2011 574

dx.doi.org/10.1021/ed100414p | J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

would still take 1014 supercomputers the age of the universe to complete. Quite a challenge! Another example is the Ising spin model, which is described fully later on. In this model, many unpaired (electron) spins interact with one another. A chemical example is the organometallic molecule FeTAC, which has long chains of Fe atoms each with unpaired electrons. In the chain of spins, each unpaired electron interacts with its nearest neighbor and the aim is to calculate the resulting energy, heat capacity, entropy, and magnetic susceptibility at any given temperature. If the number of spins is small, for instance 100, because each electron’s spin is either parallel or antiparallel to its nearest neighbor, there are 2100 or ≈1030 configurations to compute to find the energy. To calculate the partition function by direct summation is a computational disaster because most states have a similar energy, and the summation will only reach a fixed value when all terms are included. These examples illustrate that an approach that avoids calculating the partition function is required, and this is where eq 3 becomes important.

Consider next either the classical harmonic oscillator or the quantum harmonic oscillator at high temperature. There is a kinetic energy contribution to the total energy because the atoms are moving relative to one another about their equilibrium position with velocity v. Equation 7 gives the average kinetic energy. Additionally, there is the potential energy Q(x) = Ep = kx2/2, where x is the displacement from the equilibrium position. Mathematically this is the same calculation as eq 7 and produces the same result, but x is used instead of v and k instead of m.a The total average energy is therefore kBT. This result is a consequence of the equipartition theorem of classical mechanics where each quadratic term in the energy representing either the momentum px (kinetic energy E = p2x/(2m)) or position x, contributes kBT/2 to the total energy.1,2 To calculate ÆEæ, a Monte Carlo method can also be used. Just as in an experiment where repeated measurements improve the precision, in a Monte Carlo calculation increasing the number of trials acts similarly. Experimental measurements on individual molecules produce results that can be different from those on an ensemble of a large number because rare events can be observed. In an ensemble measurement, the fluctuations experienced by any individual molecule are always averaged out. These fluctuations are a result of numerous interactions with the “thermal bath”. To illustrate the difference between a single versus many measurements, the standard deviation σ of the average energy ÆEæ can be used. This is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ σ ¼ ÆE2 æ - ÆEæ2 ;

’ DIRECT CALCULATION OF THE AVERAGE ENERGY To use eq 3, the average energies ÆEæ and ÆE2æ are required. When energy is a continuous variable, its average is defined as R EPðEÞ dE ÆEæ ¼ R ð5Þ PðEÞ dE where P(E) defines the distribution P(E) = g(E)e-E/(kBT) and g(E) is the number of states with energy E to E þ dE. However, a similar relation is true for any quantity, and thus, the general expression for an average quantity of a continuous variable Q(x) is R Q ðxÞPðxÞ dx ð6Þ ÆQ æ ¼ R PðxÞ dx

and the fractional error, using eq 3 is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ kB T 2 CV ¼ ÆEæ ÆEæ

For a single particle, CV = kB/2 and ÆEæ = kBT/2, making σ/ÆEæ = 1. Thus, the average energy obtained from a measurement on a single particle has, not surprisingly, a very large standard deviation. When N particles are measured (or N repeated measurements of the same one), then CV = NkB/2, ÆEæ = NkBT/2, and now σ/ÆEæ ≈ 1(N)1/2, which is a vanishingly small quantity if N is large, and a good estimate of the true ÆEæ is obtained. Similarly, with a Monte Carlo simulation, the more samples that are included, the more accurate the result becomes and the standard √ deviation improves in proportion to 1/ N.

where P is the probability distribution appropriate to Q. The denominator is the normalizing term equivalent to dividing by N when calculating the simple average of N numbers. In calculations involving atomic and molecular properties, P is the Boltzmann distribution and the denominator is the partition function Z. If the energy levels are discrete, each integral is replaced with a summation; for example ÆQ æ ¼

∑n Qngn e-E =k T =Z n

B

To illustrate using eq 6, consider calculating the average kinetic energy of a molecule of mass m moving in one dimension, but doing so in contact with a heat bath at temperature T so that it undergoes random thermal motion. If v is the velocity at any instant, the average kinetic energy ÆEæk is found with Q(v) = mv2/2, and the ratio of integrals is R¥ 2 m -¥ v2 e-mv =ð2kB TÞ dv 1 R¥ ð7Þ ¼ kB T ÆEæk ¼ 2 -¥ e-mv2 =ð2kB TÞ dv 2 From eq 3, as ÆE2æk = 3(kBT/2)2, the heat capacity is CV = kB/2. It should be noted here that to use eq 5, the energy must always be specified in terms of its parameters, as in eq 7. In this example, this also means changing dv to dE via dv = dE/(2mE)1/2. Only then will integrating over E return the correct result.

ð9Þ

’ THE METROPOLIS ALGORITHM The Metropolis algorithm3 is a general Monte Carlo method that allows various average properties to be calculated, such as ÆEæ or ÆE2æ, by randomly sampling configurations and retaining mainly those that are important contributors to the average value.4-6 This method of sampling only important contributions is (naturally) called importance sampling.6,7 In the Metropolis Algorithm, because the Boltzmann distribution is usually used, the sampling is energy biased. To obtain the average ÆQæ, as defined in eq 6, either the two integrals can be determined separately or their ratio determined. The unique feature of the Metropolis algorithm is that it allows an estimation of the ratio of the integrals to be made in an efficient way, and hence, averages are determined directly. This avoids having to calculate Z, which, as we have seen, may be impossible. The algorithm overcomes the necessity to search all configurations by using the Boltzmann distribution of energies to 575

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

bias the random choice of which configuration to add to the total. In doing this, the algorithm tries to add to the estimate of ÆQæ only those configurations contributing significantly to its value. In effect, the algorithm performs a random walk among the configurations, which might be energy, velocity, or position. This random walk is also called a Markov process, which is defined as a process whereby a “system” goes from one state to another in a random fashion but has no memory of its previous condition. To sample points according to the Boltzmann distribution, it is sufficient, but not necessary, to impose “detailed balance” between any two configurations. This means that if w12 is the rate of going from configuration (or state) 1 with energy E1, to a new one 2 with energy E2, then the reverse transition w21 is related to it by w12 e-E2 =ðkB TÞ ¼ -E =ðk TÞ ¼ e-ΔE=ðkB TÞ w21 e 1 B

ΔE ¼ E2 - E1 ð10Þ

The ratio of rates is thus the same as the ratio of probabilities of moving between two configurations. Detailed balance, eq 10, ensures that any one configuration can be reached from any other in a finite number of steps. This behavior is called ergodicity and it means that no configurations are systematically missed. (Alternatively, this can be expressed by saying that the time average is the same as the ensemble average.) The essential step in the Metropolis algorithm is the method of choosing the ratio of probabilities. Two energies E1 and E2 from two guessed configurations, the current one 1, and a new one, 2, are calculated. If the new one has the lower energy, then this is accepted. If not, the chance of 2 contributing to ÆQæ is guessed by comparing e-(E2-E1)/(kBT) to a uniform random number in the range 0 to 1, and accepting 2 if the random number is smaller. If not, the current configuration 1 is used again.

’ OUTLINE OF THE METROPOLIS ALGORITHM The algorithm to calculate an average quantity ÆQæ is summarized as (1) Initialize parameters, for example, Qtot = 0 Calculate E1 using any reasonable initial values. Calculate the initial Q1. (2) Start a loop Calculate E2 using other parameter values chosen at random. (i) If E2 < E1, keep the new state 2 and add Q2 to a total quantity Qtot. (ii) If E2 > E1, then calculate e-(E2-E1)/(kBT): (a) If e-(E2-E1)/(kBT) > r, where r is a uniform random number between 0 and 1, keep the new state 2 and add Q2 to Qtot. (b) If e-(E2-E1)/(kBT) e r, retain the old state 1 and add Q1 to Qtot. (3) Continue the loop until Qtot is obtained to a sufficient accuracy or a predetermined number of trials are completed. (4) Average ÆQæ = Qtot/N for N calculations. Notice that even if a new state is not accepted, the last used value Q1 is still added to the total. This is because each (randomly chosen) configuration has to be added in. If a sufficient number of samples are used in the calculation, all configurations will be sampled in proportion to their importance in contributing to ÆQæ.

As an illustration, the mean displacement, Æxæ, the mean square displacement Æx2æ, the mean potential energy and the vibrational heat capacity of a classical harmonic oscillator are calculated.

’ EXAMPLE OF THE METROPOLIS ALGORITHM: HARMONIC OSCILLATOR The harmonic oscillator with force constant k has potential energy Ep = kx2/2 at displacement x from its minimum energy. The average displacement is zero, but the mean square displacement is not. Using eq 7, the average kinetic or potential energy is ÆEæ = kBT/2 and its square is ÆE2æ = 3(kBT/2)2. To make the calculation classical, many energy levels must be populated at temperature T. This implies that a very high temperature or a very small force constant should be used such that k/T is small, for example 0.05. A temperature of 300 K and k = 10 N m-1 are used. This force constant is far smaller than that for many molecules although alkali metal dimers do have small force constants. To calculate the average energy, random guesses of the values of the displacement are made. To allow complete sampling of the distribution, the maximum displacement must be large enough to make e-Ep/(kBT) small, ≈10-8 or less. In addition, the range of x must be positive as well as negative, and the random number generator must ensure this. If this is not done, then microscopic reversibility is not achieved because not all the possible configurations, in this case bond extension, can be sampled. The details of the calculation are given in the Supporting Information. A typical calculation produces 2.05  10-21 J for the average energy (Eav, see Maple code), 1.27  10-41 J2 for its square (E2av), and 6.88  10-24 J K-1 for the heat capacity (CV). These will vary slightly each time the calculation is run, but when the number of trials is large, ≈ 20,000, they are always close to the theoretical values (per molecule), which are ÆEæ = 2.07  10-21 J, ÆE2æ = 1.285  10-41 J2, and CV = 6.90  10-24 J K-1 at 300 K. The calculation of the average displacement and its square are outlined in the Supporting Information. Typical calculated values are Æxæ  Xtot/N = -0.000054 nm, which is, as expected, effectively zero, and Æx2æ  X2tot/N = 0.00209 nm2, which is close to the theoretical value of kBT/k = 0.00207 nm2. This is in agreement with the uncertainty in the positions of the atoms when a molecular structure is determined by X-ray crystallography. The atoms’ thermal motion allows Æx2æ to be measured experimentally. The thermal effect on an atom’s position is often shown in a structure as an ellipsoid surrounding each atom. As with real experiments, all these Monte Carlo calculations produce an average value ( statistical error. The size of the error is reduced in inverse proportion to the square root of the number of trials. A graph showing results of this calculation is given in the Supporting Information. The calculation to estimate the average kinetic energy of the harmonic oscillator mv2/2, is essentially the same as just described; however, the constants are interpreted in a different way. The constant k now represents the reduced mass, which is 10 amu, and x represents the velocity v in m s-1 about the origin of coordinates taken as the center of gravity of the molecule. The velocity must range from at least (4000 m s-1, if the exponential e-ΔE/(kBT) is to be ≈10-8 or less at the largest speed. A typical value of the average kinetic energy (per molecule) is found to be 2.03 ( 0.01  10-21 J, which is effectively the same result as for the previous Monte Carlo calculation of ÆEæ and again demonstrates the equipartition theorem.1,2 576

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

this interaction to the temperature. A small section of the spin state is shown in Figure 1B. Parallel spin pairs have energy -J but þJ if they are antiparallel. If J is positive, the spin system is ferromagnetic and if negative, anti-ferromagnetic. The energy of a chain in the absence of an external magnetic field is given by ER ¼ - J

∑k sksk þ 1

ð12Þ

where sk = (1 represents the spin on the kth site, and the sum is taken over all pairs skskþ1. The various combinations the spins take produce the different energies ER. Because the exchange interaction depends on the overlap of wavefunctions, as does the related Coulomb interaction, it is of very short-range; therefore, only the nearest neighbors are important, and the summation of eq 12 is pairwise. Before doing the calculation, it is worth discussing with students in general terms how ÆEæ, CV, and S are expected to vary with temperature as this may provide some guide to their own calculations. Because the spins interact with one another, at very low temperature they will all align in one direction with a total energy of -NJ (assuming J > 0) as there is not enough thermal energy to disrupt spin alignment. As the temperature increases, spins will flip over. Eventually, at infinite temperature, which in practice may be well below room temperature, the energy approaches zero because there are as many parallel as antiparallel spin pairs. The spin heat capacity at low temperatures is also close to zero, because there is not enough energy to flip a spin, and therefore, the rate of change of energy with temperature is zero; limTf0CV(T) = limTf0(∂U/∂T)V = 0. At infinite temperature, the heat capacity again approaches zero. Because there are now equal populations of parallel and antiparallel spin pairs, no more energy can be absorbed, and thus, the rate of change of energy with temperature is zero. As CV is zero at both extremes of temperature and is, by definition, a positive quantity, it must pass through a maximum at some finite temperature whose value will depend on J. The entropy of the chain of N spins is kBln(2) at zero temperature, because two degenerate and ordered spin states exist. At high temperatures, where the thermal energy is far greater than the spin interaction energy, the entropy is NkB ln(2). This is N times greater than at 0 K and a dependence on N is to be expected because entropy is an extensive quantity. The entropy therefore initially increases with an increase in temperature and then becomes constant. The magnetization and susceptibility can also be simulated. Details of these calculations are given in the Supporting Information. Without having to evaluate the partition function directly, the energy, heat capacity, and entropy can be calculated using the Metropolis algorithm. In the calculation shown in Figure 2, a chain of 500 spins is simulated and the number of spin states is therefore 2500 ≈ 3  10150, all of which would have to be evaluated in a direct calculation. The results of the Monte Carlo calculation, Figure 2, have been obtained by randomly sampling each site 4000 times. Although this represents a large number of calculations, it is only a miniscule fraction of the total number of possible spin states and demonstrates the power of this algorithm. (The code for the algorithm is in the Supporting Information.) The solid lines, Figure 2, are calculated with the analytical solution to the partition function whose derivation is also outlined in the Supporting Information. It is clear that the

Figure 1. (A) Schematic of part of the FeTAC crystal structure, based on Greeney et al.8 The Fe atoms (center) are linked by two chlorines (gray), and the oxygen atoms (black) complete the coordination. In this compound, J/kB = 8.85 K. (Note that Landee et al.9 define the energy with 2J instead of J as shown in eq 12 and hence report J/kB = 17.7 K.) (B) Representation of electron spins as either “up” or “down”.

If the kinetic energy mv2/2 is interpreted as that due to molecules of a gas colliding with one another rather than that of the atoms of a harmonic oscillator, then the Maxwell speed distribution can be constructed. To do this, the velocities that are accepted in the Metropolis step are stored as the calculation proceeds. Velocity is a vector but speed is not; consequently, the angles that a molecule has with respect to any axes are integrated away on calculating the speed distribution and no longer enter into the calculation. The result is that the speed distribution depends only on v2, which is interpreted as a scalar equal to the speed squared. For example, the probability per unit speed of molecules with speed s = |v| is,   4 m 3=2 2 -ms2 =ð2kB TÞ se ð11Þ f ðsÞ ¼ pffiffiffi π 2kB T R It can be verified that ¥ 0 f(s) ds = 1. If the values of v calculated from the Monte Carlo method are put into eq 11, they are distributed toward smaller values of v, which shows the weighting that the method produces. However, if the Metropolis probability test is changed to (s22/s21)e-ΔE/(kBT) from e-ΔE/(kBT), then the results are distributed according to the Maxwell distribution. A plot is shown in the Supporting Information.

’ ONE-DIMENSIONAL SPINS AND THE ISING MODEL In the calculations described so far, no interactions occur between the particles. One of the simplest examples in which interaction is allowed is the one-dimensional (1-D) Ising spin model. This is realized in practice with some Fe(II)Cl3 compounds which, in the crystalline state, have chains of Fe2þ atoms held apart with dichloride bridges.8,9 Part of the crystal structure of trimethylammonium iron(II) trichloride hydrate, (FeTAC), [(CH3)3NH]FeCl3 3 2H2O, is shown in Figure 1. The coupling between different chains of molecules in the crystal is relatively weak, but the interaction along a chain is so strong that this compound behaves as a 1-D chain of interacting electron spins. In the 1-D Ising model, interacting spins are aligned in a chain. A spin may flip because of its interaction with surrounding atoms and consequently neighboring spins are affected. These spins in turn affect their neighbors, and thus, one flip may propagate along the chain until all spins come to equilibrium with their surroundings. The electrons interact with one another by the exchange interaction of magnitude J. The number of electrons in each spin state, either “up” or “down”, depends on the ratio of 577

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

Figure 4. (A) Examples of spin states for 500 spins and 2000 samples/ spin and with J/kB = 1 K at T = 1/2, 1, and 10 K. The pattern at T = 1/2 K contains 11 groups since the black areas correspond to one spin state and gray the other. (B) The open circles show the Monte Carlo calculated spin entropy as ÆSæ/kBN vs temperature with J/kB = 1 K, J/kB = 5 K and using eq 13. At each temperature and J value, N = 500 and 2000 samples/spin were taken. The solid lines are calculated with the theoretical function (see the Supporting Information, eq S7). The high temperature limit for both J values was calculated at T = 100 and was effectively kBNln(2). This limit is indicated by the dashed line which has a value of ln(2).

calculated outside the Monte Carlo loop, which saves considerable time. The second consideration concerns the finite length of the chain. Theoretically, infinite length is assumed, but obviously this must be finite and thus a boundary condition is necessary to simulate infinity. The boundary condition used effectively makes the chain circular, so that if a spin’s nearest neighbor is one step off the end of the chain, the state of the spin at the other end is used. Therefore, if the indices run from 1 to N, then index N þ 1 f 1 and index 0 f N. These steps are labeled in the Maple code shown in the Supporting Information. The last consideration concerns the initial arrangement of the spins. Are these to be aligned, to be random, or somewhere in between? Although it does not matter how the initial spins are arranged provided enough Monte Carlo runs are taken, at low temperatures fewer runs are needed if the spins are initially parallel.

Figure 2. The open circles are (A) the Monte Carlo calculated energy, ÆEæ, per spin; ÆEæ/JN. (B) The heat capacity per spin ÆCVæ/kBN vs temperature (K), with J/kB = 1 K. One run of 4000 samples is shown for each point. The average of five runs of ÆCVæ/kBN at T = 1 K has a mean of 0.45 and a 95% confidence limit of (0.015. The lines are from the analytical functions; see the Supporting Information.

Figure 3. A schematic of two configurations of spins related by a single spin flip. Only values of ΔE = 0, (4J occur. The interaction energy is -J between a pair of parallel spins whether up or down, and þJ if antiparallel. The energy change is calculated as final - initial state.

Metropolis algorithm accurately simulates the change in energy and heat capacity with temperature. In implementing this particular Monte Carlo calculation, there are three considerations. The first is to calculate the energy change caused by flipping one spin, the second is to decide what happens at the ends of the chain, and the third is to determine what values the spins should have at the start of the calculation. First, to determine the energy change, a site is chosen at random and its spin flipped. The change in energy this causes is found by looking at the site and its two neighbors and calculating what happens; Figure 3 shows two examples. The change in energy ΔE = E2-E1 is either 0 or (4J, and because a negative or zero change in energy is always accepted, the exponential expression in the Metropolis sampling step is always e-4J/(kBT). This can be

’ CALCULATING THE ENTROPY Although the entropy at a given temperature T can be calculated by integrating CV/T, it can also be calculated from the arrangement of the spins. Maps of the spin state at the end of the calculation indicate how the entropy changes with temperature, and three examples are shown in Figure 4A. Informally, we may label a black group as “parallel spin up” and a white one “parallel spin down”, although “up” and “down” should not be taken literally. Within any blocks of either color, the spins are parallel to one another; hence, the coupling is -J between adjacent spins. A color boundary thus indicates spin pairing and a coupling of þJ between these two spins; see Figure 3. At T = 0 K, all spins are in the same state, which would mean that the whole strip would be of the same color. However, there are two ground states at T = 0 and these could be colored either black or white. At low temperatures, for 578

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

example, T = 1/2 or 1 K, (Figure 4A) spins are grouped into a few large blocks of similar spin state. Although only one example is given in the figure, it is clear that if many calculations are performed at, say, T = 1/2 K, only a limited number of arrangements of these large blocks of spins will be possible. These configurations or arrangements determine the entropy. As the temperature is increased, the coupling between spins becomes relatively less important compared to the energy supplied by the surrounding heat bath. This reduces the size of the groups of similar spin state but increases their number, which leads to more of ways of arranging them, and hence, a larger entropy. The entropy of the equilibrium state is given by S = kB ln(Ω) where Ω is the number of arrangements that the spins can take. To calculate Ω, the number of groups of spins p present at the end of the calculation and at a given temperature is counted; see Figure 4A. This is done by searching the pattern of spins for the number of spin changes and adding one to the total (see the Supporting Information, eq S7). From the theory of statistics, the number of ways of choosing Np distinguishable objects from a group of N distinguishable objects, where the order of choosing does not matter, is given by Ω0 ¼

N! Np !ðN - Np Þ!

ð13Þ

This is the number of arrangements, or combinations. Equation 13 is a good approximation to the number of arrangements of groups of spins, when N is large and the temperature is not zero, thus Ω ≈ Ω0 . The good agreement, Figure 4B, between the equation describing the entropy per spin (see the Supporting Information) and that obtained using eq 13 shows that this is a good approximation. When T f ¥, the number of spin arrangements must be greatest because the interaction between spins is unimportant compared to thermal energy. From eq 13, this occurs when Np = N/2, and provided N is large Ωmax f 2N. The maximum entropy is therefore Smax = NkBln(2) and the entropy per spin when T f ¥ is kB ln(2). This is the same result as obtained from theory (see the Supporting Information). Because the entropy depends on N, this result demonstrates that entropy is an extensive quantity. At T = 0, all the couplings between spins will be the same (at -J) and one group of spins exists, and thus, a value of Ω = 1 might seem appropriate and this would make the entropy zero. However, because of the quantum properties of electron spin, two ground states exit and because both of these are populated at T = 0, Ω = 2 and the total entropy is kBln(2). When N f ¥, the entropy per spin tends to zero, S0 = kB ln(2)/ N f 0, in accordance with the third law of thermodynamics. Within statistical variation, the entropy per spin obtained by counting groups of spins and using eq 13 and that calculated from theory (Supporting Information) are the same, see Figure 4B. Five repeated measurements at J/kB = 1 K produced ÆSæ/kBN = 0.34 ( 0.04 (95% confidence limits) at T = 1 K, and of 0.58 ( 0.03 at T = 2 K. Finally, we note that S = kB ln(Ω) is derived using a microcanonical ensemble, (N, V, E constant), whereas the entropy calculated using the partition function is derived from a canonical ensemble (N, V, T constant). Energy fluctuations are permitted in the canonical ensemble but energy is fixed in the microcanonical ensemble. This inherent difference does not contradict their equivalency, provided that the size of the fluctuations becomes vanishingly small in the canonical ensemble, which it does in the

limit of large N. The argument to demonstrate this follows along the lines of that of eq 9 where the energy dispersion σ (eq 8) is now interpreted as the energy fluctuation of the canonical ensemble. The ratio of energy fluctuation to average energy σ/ÆEæ ≈ 1/(N)1/2, becomes vanishingly small as N increases.10

’ PHASE CHANGE The gradual and continuous change in the heat capacity and entropy indicates that the phase transition between ordered and disordered spins is unlike the familiar first-order or “all-or-none” phase transitions1,2 or even that in the two-dimensional Ising spin model.6,7 Whereas the general change in CV and S with temperature is similar to that of a normal substance (O2, H2O, etc.), the detail is quite different. This is because in the 1-D spin case the heat capacity increases and decreases in a continuous manner with increasing temperature. Similarly, the entropy smoothly increases with temperature to reach its final value. In a normal substance, however, there is a step change both in CV and S at a phase boundary.1,2 The cause of the different behavior is because, energetically, it is more favorable for two regions of opposite spins to co-exist than it is for either one or the other to exist on its own. Interestingly, a similar behavior occurs with the co-existence of helix and coil regions in proteins.11 ’ CONCLUSIONS A Monte Carlo method, called the Metropolis algorithm, has been described and its use illustrated as being effective in calculating thermodynamic quantities such as energy, heat capacity, and entropy. Results for two different one-dimensional problems, the harmonic oscillator, and the Ising spin model are presented. The Monte Carlo calculations match exactly the analytic functions within statistical limits with no free parameters. Although the method is generally useful, it is particularly so when the partition function contains far too many terms to be calculated directly. The method illustrated should enable undergraduates to gain an insight into the statistical nature of thermodynamic properties as well as an alternative way of numerically calculating quantities. Annotated code based on the symbolic-algebra language Maple is given in the Supporting Information. The programs are written using expressions common to most computer languages, so that it should be straightforward to convert them, for example, for use in Mathematica, MathCAD, MATLAB, Basic, or Fortran 77. ’ ASSOCIATED CONTENT

bS

Supporting Information Maple Code for the Metropolis algorithm, The Harmonic Oscillator, Average Displacement, Maxwell Speed Distribution 1-D Ising Spins, Entropy, Magnetization and Susceptibility. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT I am grateful to G. Reid and M. deMiranda, (School of Chemistry, University of Leeds) for several helpful discussions. 579

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580

Journal of Chemical Education

ARTICLE

’ ADDITIONAL NOTE a Although x and v are dummy variables in the integration, it is useful to label them separately. ’ REFERENCES (1) Atkins, P.; de Paula, J. Physical Chemistry; 8th ed.; Oxford University Press: Oxford, 2006. (2) Engel, T.; Reid, P. Physical Chemistry; Pearson: San Francisco, CA, 2006. (3) Metropolis, N.; Rosenbluth, A. M.; Teller, A.; Teller, E. J. Chem. Phys. 1953, 21, 1087–1092. (4) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 1996. (5) Newman, E.; Barkema, G. Monte Carlo Methods in Statistical Physics; Clarendon Press: Oxford, 1999. (6) Krauth, W. Algorithms & Computations; Oxford University Press: Oxford, 2006. (7) Beddard, G. S. Applying Maths in the Chemical & Bio-molecular Sciences; Oxford University Press: Oxford, 2009. (8) Greeney, R. E.; Landee, C. P.; Zhang, J. H.; Reiff, W. M. Phys. Rev. B 1989, 39, 12200–12214. (9) Landee, C. P.; Kuentzler, R.; Williams, J. J. M. J. Appl. Phys. 1990, 67, 5604–5606. (10) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: Oxford, 1987. (11) Finkelstein, A. V.; Ptitsyn, O. B. Protein Physics, A Course of Lectures.; Academic Press: San Diego, CA, 2002.

580

dx.doi.org/10.1021/ed100414p |J. Chem. Educ. 2011, 88, 574–580