MIT 3.320/SMA 5.107/CMI Atomistic Modeling of Materials Spring 2005Lab 5: Monte Carlo simulationsDue date: 5/05/2005In this lab, we will be using a Monte Carlo code that we have written ourselves. MonteCarlo codes are usually simple enough so that you can write them yourself.The problem: Adsorption of H on the (001) Surface of PdWhen hydrogen adsorbs onto a clean (001) surface of Pd, the H atoms sit between the Pdatoms. These possible adsorption sites are marked with an X in Figure 1. The Hcoverage, θ, is the fraction of possible adsorption sites that are filled. Evidence obtainedwith LEED [Behm et al, Surface Science 99 (1980) p320] shows that for θ = 0.5 the Hatoms form an ordered arrangement. We will model this adsorption problem with asquare 2D lattice model that consists of all the sites where H can adsorb.Figure 1. Lattice model for H adsorption on Pd.The energy modelAs you learned in class, Monte Carlo trajectories are determined by the Metropolisalgorithm. If ∆E<0 for a trial perturbation, the simulation accepts the perturbation; if∆E>0 for a trial perturbation, the simulation accepts the perturbation with probability P α exp(-∆E/kT). All Monte Carlo simulations need an energy model to calculate ∆E. Inprinciple, if you really wanted to, you could use a first-principles calculation to calculate∆E for every trial perturbation. In practice, this takes way too long, and simpler energymodels are required. The energy model we will be using is called the "Ising Model". The figure belowshows a two-dimensional square lattice on which a nearest-neighbor (NN) and next-nearest neighbor (NNN) interaction are defined. The occupation on each lattice site i ischaracterized by a spin occupation variable, σi, which can take the value -1 or +1. In amodel for a binary alloy, +1 and -1 can respectively stand for occupation of that latticesite by an A or B atom. In this problem, σi =+1 if the site is has an H adsorbed; σi = -1 ifthe site is vacant.1XXXXXXXXXXXXP d a to m sH a d s o rp tio n s ite sMIT 3.320/SMA 5.107/CMI Atomistic Modeling of Materials Spring 2005Figure 2. Definition ofinteractions characterized byV1 and interactionscharacterized by V2.The Hamiltonian (energy expression) of this lattice model is given below:H =12V1σii, jNN∑σj + 12V2σii, jNNN∑σj − µ σiiN∑The first summation is performed over first-nearest neighbor pairs; the second summationis performed over second-nearest neighbor pairs. The factor of 1/2 accounts for the factthat every pair is counted twice in the sum. V1 and V2 are the first and second neighborinteractions. µ is the chemical potential (or magnetic field): it fixes the amount of Hadsorbed. At a given temperature it is linearly related to the log(PH2). Note that θ = (1 + <σ> )/2, where <σ> is the average of the spin occupation variable.If the parameters V1 and V2 are chosen well, and enough neighbors are included inthe model, the Ising model can do an excellent job of representing the actual energetics ofthe system. A first-principles tight binding study found that V1 is positive and V2 isnegative, with |V2| > V1. We will therefore use the ratio V2/V1 = -2 to model thissystem. Grand canonical simulationsMonte Carlo simulations on lattice models are typically performed in the grand-canonicalensemble (Fixed volume, temperature, and chemical potential - composition canfluctuate) rather than the canonical ensemble (fixed volume, temperature, andcomposition - chemical potential can fluctuate). This is largely for computationalefficiency. In a canonical simulation, picking a trial perturbation is slow: for every trialspin you flip, you have to flip another one (after making sure it is the correct spin) topreserve concentration. Both canonical and grand canonical ensembles can be shown togive the same result for thermodynamic averages. As you may be used to thinking of composition as an independent variable, thismay require you to change your thinking of how you plot thermodynamic quantities. In2V1V2MIT 3.320/SMA 5.107/CMI Atomistic Modeling of Materials Spring 2005the Monte Carlo simulations we do not pick the concentration directly: we pick µ, and thesystem will find an equilibrium concentration. Also, using the grand canonical ensemble changes the Metropolis algorithmsomewhat. Instead of looking at ∆E, we use ∆E-µσi, where σi is the site of the trialperturbation.Running the codeThe executable code is contained within a java archive file (2DMC-1.0.jar). To run thecode you will need to start up the java virtual machine (JVM) with the jar file in theJVM's classpath and call the RunMC program, i.e.user@hpcbeo2: java -cp 2DMC-1.0.jar RunMC “input_filename”The input file contains tokens and values (one per line). To see a description of theavailable (and necessary) tokens and values execute the command:user@hpcbeo2: java -cp 2DMC-1.0.jar RunMC -hMore Helpful hints:Before you start calculations it is important to get a sense of how long it will take you toobtain a certain result. This may be part of the input you need to decide whether youactually should do the calculations ! Never start up a calculation unless you have someidea whether it will take minutes, hours or days to finish.µ-range to consider. To get an idea of the chemical potentials in your system, you caneither do a few quick short runs, or (better) look at the energy of some possible structuresin your system and see where they would cross. Τ-range to consider. To get an idea of where the transition(s) is (are), you should do afew quick runs over a coarse temperature range. Once you determine the range ofinterest, you can do runs with a finer temperature mesh.How do I look for the phase transitions?A simple way you can look for first-order transitions is to plot the chemical potential vs.spin. Run the Monte Carlo simulation at a fixed temperature, but vary ?. Fromthermodynamics, we know that in a 2-phase region, the chemical potential will beconstant. Thus, at a fixed temperature, if we plot chemical potential vs. spin, we can findthe boundaries of the two-phase region. Be careful! This applies to first-ordertransitions, but you should not rely upon this method to find boundaries of second-ordertransitions.3MIT 3.320/SMA 5.107/CMI Atomistic Modeling of Materials Spring 2005In lab 4, you learned that one can find a first-order transition by looking