MacKinnon - Computational Physics, страница 11

These are useful for describing isothermal or isobaric processes.Note that for a set of mutually attractive particles it may not be necessary to constrain the volume butfor mutually repulsive particles it certainly is necessary.4.6 Problemsœ ÷  n1. Show that (4.1) has the property that it generates all the integers from toorder ifand, and that the sequence repeats itself thereafter.in an apparently random2. A certain statistical quantity is distributed according to\"<Ky&; f K"]KÑ&for —lK — Given a function which generates random numbers uniformly distributed betweento transform these into the distribution \. & , show howThe Ising Model43S> e ¿ e ÿ ª ÿ e ª 1K3. Suggest a Monte–Carlo integration procedure for the integral_`†…„‡Rˆ ˆ„`4.

Describe how you would use the Metropolis algorithm to generate a set of random numbers distributed according to the integrand of problem 3.The following are essay type questions which are provided as examples of the sorts of questions whichmight arise in an examination of Monte–Carlo and related methods.5. In an ionic conductor, such as AgI, there are several places available on the lattice for each Agion, and the ions can move relatively easily between these sites, subject to the Coulomb repulsionbetween the ions.Describe how you would use the Metropolis algorithm to simulate such a system." & ¿e6. Some quantum mechanics textbooks suggest that the ground state wave function for a Hydrogenatom is37 ;‰r[Describe how you would use the variational quantum Monte–Carlo procedure to calculate the groundstate energy for a range of values of and hence how you would provide estimates of the true groundstate energy and the value of .7.

Describe how you would simulate the melting (or sublimation) of Argon in vacuo using the moleculardynamics method. Pay particular attention to how you would calculate the specific and latent heats.4.7 Project — The Ising Model4.7.1 IntroductionÎ €  €€šºL$ …΀The Ising model for a ferromagnet is not only a very simple model which has a phase transition, but it canalso be used to describe phase transitions in a whole range of other physical systems. The model is definedusing the equationsS S(4.27)Òâthe valueswhere the and designate points on a lattice and S takesssystems differ in the definition and sign of the various’s.T( 0 .

The various different physical4.7.2 The Model and Method €€º$ … á š ÞÞ Òwhere â is only summed over the 4 nearest neighbours of .This model can be studied using the Metropolis method as described in the notes, where the state canbe changed by flipping a single spin. Note that the change in€ energydue to flipping the Ë th spin from Š toŒ‹is given byŒ€,º L$ s … S [(4.29)ÍHere we will consider the simple case of a 2 dimensional square lattice with interactions only betweennearest neighbours. In this casesS S(4.28)The only quantity which actually occurs in the calculation isSŽn QÍhJ!w"%$ ,º 7 Ë & :Í(4.30)The Ising Model44‹and this can only take one of five different values given by the number of neighbouring spins.

Hence it issensible to store these ins a short array before starting the calculation. Note also that there is really only 1n Q , so that it would make sense to write your program in terms of this singleparameter in the model,sparameter rather than and Q separately.The calculation should use periodic boundary conditions, in order to avoid spurious effects due toboundaries. There are several different ways to achieve this. One of the most efficient is to think of thesystem as a single line of spins wrapped round a torus. This way it is possible to avoid a lot of checking forthe boundary.

For an system of spins define an arrayelements using the shortest sensible‹ of variable type: char in C(++). Itiseasiertouseforspinandforspin Š , as this makes the calculation‹of the number of neighbouring spins easier. In order to map between spins in a 2d space r and in the 1dthe following mapping can be used.arrayy7 Ë¥ßf 0ÍßS F› S S#T : Sr e x F› S S ÿ egT : SrS y F› S S : Sr e y F› S S ÿ e (4.31)ë0 Íë Íë Íë 0Íwhere the 2nd elements of the array are always maintained equal to the 1st . This way it is nevernecessary to check whether one of the neighbours is over the edge. It is important to remember to changeSS ÿ whenever S is changed.ÍÍSr xThe calculation proceeds as follows:1. Initialise the spins, either randomly or aligned.2. Choose a spin to flip. It is better to choose a spin at random rather than systematically as systematicchoices can lead to spurious temperature gradients across the system.3.

Decide whether to flip the spin by using the Metropolis condition (see notes).4. If the spin is to be flipped, do so but remember to flip its mirror in the array.5. Update the energy and magnetisation.6. Add the contributions to the required averages.7.

Return to step 2 and repeat.4.7.3 The PhysicsIn general it is advisable to run the program for some time to allow it to reach equilibrium before tryingto calculate any averages. Close to a phase transition it is often necessary to run for much longer to reachequilibrium. The behaviour of the total energy during the run is usually a good guide to whether equilibriumhas been reached. The total energy, , and the magnetisation can be calculated from (4.27) andº… ßa ¥ (4.32)It should be possible to calculate these as you go along, by accumulating the changes rather than by recalculating the complete sum after each step.

Alattice should suffice for most purposes and certainlyfor testing, but you may require a much bigger lattice close to a transition.A useful trick is to use the final state at one temperature as the initial state for the next slightly differenttemperature. That way the system won’t need so long to reach equilibrium.It should be possible to calculate the specific heat and the magnetic susceptibility. The specific heatcould be calculated by differentiating the energy with respect to temperature. This is a numerically questionable procedure however.

Much better is to use the relationshipes ‰‘fj(4.33)in Q0 º0 $ º 0ŒË  °±Quantum Monte Carlo Calculation45eSimilarly, in the paramagnetic state, the susceptibility can be calculatedusingsvfj’ ‰in Qß 0 $ ß 0±°(4.34) ËßÊ ä ß and the averages are over different states, i.e. can be calculated by averaging over thewheredifferent Metropolis steps. Both these quantities are expected to diverge at the transition, but the divergencewill tend to be rounded off due to the small size of the system. Note however that the fact that (4.33) &(4.34) have the form of variances, and that these diverge at the transition, indicates that the average energyand magnetisation will be subject to large fluctuations around the transition.Finally a warning.

A common error made in such calculations is to add a contribution to the averagesonly when a spin is flipped. In fact this is wrong as the fact that it isn’t flipped means that the original statehas a higher probability of occupation.4.8 Project — Quantum Monte Carlo Calculation4.8.1 IntroductionThis project is to use the variational quantum Monte Carlo method to calculate the ground state energy ofthe He atom. The He atom is a two electron problem which cannot be solved analytically and so numericalmethods are necessary. Quantum Monte Carlo is one of the more interesting of the possible approaches(although there are better methods for this particular problem).4.8.2 The MethodThe Schrödinger equation for the He atom in atomic units is,f f Œ $ f Ý 0; Z $ f Ý 0; ÿ $ r T $ r 0 r T 0  w " rT : r0 &;º w " rT : r0 & :(4.35)TT z rT z , r 0 z r0 z , and r T 0 z rTÁ$ r0 z .where r and r 0 are the position vectors of the two electrons, rEnergies are in units of the Hartree energy (1 Hartree = 2 Rydbergs) and distances are in units of theBohr radius.

The ground state spatial wavefunction is symmetric under exchange of the two electrons (therequired antisymmetry is taken care of by the spin part of the wavefunction, which we can forget aboutotherwise).z ƒr r , is,The expression for the energy expectation value of a particular trial wavefunction,z ƒ x z ƒƒ! 9rrz ƒ z ƒ(4.36)!rrº ¬û [ [[[[ û[û û‡ T^‡ 0 [ ‡ TJ ‡ 0" T: 0&In the variational Monte Carlo method, this equation is rewritten in the form,º ƒ  >[[ [^> Œ z ƒwhere oj " rT : r0 & ‡ r T ‡ r 0 :jk" rT : r0 i& z [ƒ! [ [ " rT z : rƒ 0 z& z ƒ ƒ " ‡ rT T : r0‡ &û û ! r r09x z ƒ(4.37)(4.38)is interpreted as a probability density which is sampled using the Metropolis algorithm.

Note that theMetropolis algorithm only needs to know ratios of the probability density at different points, and so thenormalisation integral,z ƒ z ƒ!(4.39)rr> [ [[ >always cancels out and does not need to be evaluated.‡ T‡ 0 :Quantum Monte Carlo Calculation46The mean of the values of the “local energy”,x z ƒz ƒ 9:(4.40)at the various points along the Monte Carlo random walk then gives an estimate of the energy expectationvalue. By the variational principle, the exact energy expectation value is always greater than or equal tothe true ground state energy; but the Monte Carlo estimate has statistical errors and may lie below the trueground state energy if these are large enough.