Lattice-switch Monte Carlo: the fcc—bcc problem

Lattice-switch Monte Carlo is an efficient method for calculating the free energy difference between two solid phases, or a solid and a fluid phase. Here, we provide a brief introduction to the method, and list its applications since its inception. We then describe a lattice switch for the fcc and bcc phases based on the Bain orientation relationship. Finally, we present preliminary results regarding our application of the method to the fcc and bcc phases in the Lennard-Jones system. Our initial calculations reveal that the bcc phase is unstable, quickly degenerating into some as yet undetermined metastable solid phase. This renders conventional lattice-switch Monte Carlo intractable for this phase. Possible solutions to this problem are discussed.


Introduction
A common problem which occurs in condensed matter physics is as follows: for a given substance, which of two candidate phases is preferred at a given, say, temperature T and pressure P ? This problem amounts to evaluating the free energy difference between the two phases; the preferred phase has the lower free energy. Unfortunately, calculating the free energy difference between two phases to a sufficient accuracy to solve this problem can be difficult. This is the case for a plethora of systems of practical interest, and is by no means limited to 'realistic' models of particle interactions. For instance, while the hard-sphere solid is an archetype of a 'simple' system, until relatively recently there was contention regarding whether the equilibrium phase was fcc or hcp [1]. A similar situation also existed for the Lennard-Jones solid at low temperatures and pressures [2].
The method which resolved the aforementioned hard-sphere and Lennard-Jones disputes is lattice-switch Monte Carlo (LSMC) [1,2]. LSMC allows the free energy difference between two phases to be calculated efficiently. Furthermore it is 'exact' in the sense that it relies upon no approximations other than those present in the model of particle interactions it is used in conjunction with. LSMC has been applied to a wide range of systems since its inceptionas summarised in Table 1 1 . This reflects the generality of the method; it can in principle be applied to any pair of phases, and any model of particle interactions. This feature of LSMC, in combination with its supposedly superior computational efficiency compared to alternative methods, makes it an attractive prospect for ab initio applications. However, it should be noted that claims of its superiority have proved contentious, and remain somewhat of an open question [4,5]. The layout of this work is as follows. In the following section we provide a brief introduction to LSMC. Detailed accounts of the method can be found in the references. In Sec. 3 we consider the application of LSMC to the fcc and bcc phases of the Lennard-Jones solid. We first describe a lattice switch between these phases, and conclude by presenting the results of our initial investigations.

Lattice-switch Monte Carlo: an introduction
Consider the Gibbs free energy difference ∆F ≡ F 1 − F 2 between two crystalline phases 1 and 2, where F α denotes the free energy of phase α. It can be shown that where β denotes the thermodynamic beta, and p α denotes the probability of the system being in phase α at thermodynamic equilibrium -assuming that the system is constrained to be in either phase 1 or phase 2. The above equation can be exploited to calculate ∆F theoretically: extract p 2 /p 1 , the time the system spends in phase 2 relative to phase 1, from an N P T simulation of the system, and substitute this quantity into the above equation. However, this approach is intractable for simulations utilising 'realistic' particle dynamics if the typical time taken for the system to transition between the two phases is very long, in which case a reasonable estimate of p 2 /p 1 cannot be deduced in a reasonable simulation time. It may even be the case that, regardless of the phase in which the simulation is initialised, the system never transitions to the 'other' phase during the course of the simulation. This stems from the fact that, at equilibrium, the most probable states comprise two 'islands of stability' in phase space: one within phase 1 and the other within phase 2. However, these two islands are separated by an entropic barrier : a region of phase space comprised of states which are very improbable at equilibrium. Hence to transition between the phases the system must traverse the entropic barrier, the success of which is very unlikely, and hence a rare occurrence. The traditional implementation of Metropolis Monte Carlo [17] belongs to the aforementioned class of simulations which utilise 'realistic' particle dynamics. In this approach (for an N P T simulation) at each time-step a trial state of the system is generated from the current state by either altering the position of one of the particles, i.e., a 'particle move' is attempted, or altering the dimensions of the simulation box, i.e., a 'volume move' is attempted. The trial state is then accepted or rejected as the new state of the system for the next time-step according to the Metropolis algorithm [17]. The end result is that for a long simulation the states are sampled from the probability distribution corresponding to thermodynamic equilibrium. However, the important properties of the Metropolis algorithm do not rely upon the mechanism to generate trial states just described; one has considerable freedom with regards to how trial states are generated. The prospect therefore exists of generating states in a manner such that the system traverses a path in phase space which allows ∆F to be accurately calculated in relatively few time-steps. Such a path would involve frequent transitions between both phases by 'jumping over' the entropic barrier. This is what is done in LSMC; a new type of move, a lattice switch, is introduced in order to supplement the aforementioned particle and volume moves. Intuitively, the system is in phase α if the positions of the particles approximately form the crystal lattice characteristic of α at the current volume. With this in mind, if the system is in phase α, one can express the position r i of particle i as where R α,i denotes the position of the α crystal lattice site which is closest to i, and u i denotes the displacement of i from that lattice site. The trial state σ ′ generated by a lattice switch from a state σ in phase 1 shares the same set of particle displacements {u i } as σ, but the underlying set of lattice vectors for In other words, the lattice switch amounts to the following transformation for all i: Hence every time a lattice switch is accepted the system transitions directly to the 'other' phase, bypassing the entropic barrier. This is, however, only half of the story. One might expect that by regularly attempting lattice switches the system will frequently transition between phases, allowing ∆F to be efficiently evaluated as described above. However, if one does this with Metropolis Monte Carlo, one finds that lattice switches are too rarely accepted for this approach to be useful. The reason for this is that, for the states visited during a typical simulation, σ ′ is almost always of a much higher energy than σ, and hence lattice switches will almost always be rejected by the Metropolis algorithm. Crucially, there exist states for which σ and σ ′ are of comparable energy; from such states switches have a reasonable chance of success. We refer to these states as gateway states, since they provide the key to jumping between both phases. Unfortunately, these states are almost never visited dduring a Metropolis Monte Carlo simulation. To encourage successful lattice switches, we therefore use multicanonical Monte Carlo [18,19,20] -which provides a means of sampling selected states more (or less) frequently than is the case at equilibrium, while still allowing equilibrium properties of the system to be evaluated -to more frequently visit gateway states. The result is that lattice switches are accepted reasonably often, and both phases are explored in a reasonable simulation time. To elaborate, we introduce a quantity M which characterises how 'gateway-like' a state is, with M = 0 corresponding to 'perfectly gateway-like', and |M| ≫ 0 corresponding to 'very un-gateway-like'. The specific definition of the quantity M depends on the system under consideration. For solid-solid free energy differences in soft-potential systems, the following definition of M has been used: where E({r i }) denotes the energy of the state with particle positions {r i }. The first term on the right-hand side is the energy associated with the displacements {u i } for phase 1, and the second term is the analogous quantity for phase 2. Note that M({u i }) = 0 if the energies associated with {u i } in both phases are identical. In this case the energy cost of a lattice switch from either phase is zero, and hence the states associated with {u i } in both phases are gateway states. Defining a macrostate as a collection of all states with the same M, we then sample all macrostates in our multicanonical simulation with equal probability. The result is that the gateway-like macrostates are frequently visited.

The fcc-bcc transition in the Lennard-Jones solid
We now turn to the problem of using LSMC to evaluate the free energy difference between the fcc and bcc phases in the Lennard-Jones system. Our motivation behind this is twofold. Firstly, the fcc-bcc transition is of profound importance to metallurgy. Our ultimate aim is to apply LSMC to this transition using more realistic models of metals than the Lennard-Jones model, such as the embedded atom model [21], or even ab initio models. However, given that there has yet to be a LSMC study of the fcc-bcc transition, it is sensible to 'tread carefully' and first study the fcc-bcc transition using the simpler, and better understood, Lennard-Jones model before proceeding to more uncharted waters. Secondly, there has been speculation that there is a region in the phase diagram of the Lennard-Jones system at high temperatures and pressures, below the melting curve, where the bcc phase will be preferred over the fcc [22]. It would be interesting to test this hypothesis, the confirmation of which would have far-reaching consequences given the widespread use of the Lennard-Jones model to describe real systems.

The fcc-bcc lattice switch
After deciding to apply LSMC to a certain system, the first problem one encounters is the choice of lattice switch. A lattice switch is a one-to-one mapping of particle positions in one phase to another phase. Hence it is necessary that the supercells used to represent both phases have the same number of particles. The Bain orientation relationship (see, e.g., Ref. [23]) provides a means for achieving this for the case of the fcc and bcc phases. To elaborate, both fcc and bcc crystals can be recast as bct crystals: the bcc crystal is equivalent to a bct crystal in which the bct unit cell has equal dimensions, i.e. a = b = c; the fcc crystal is equivalent to a bct crystal in which a = b = c/ √ 2. This is illustrated in Fig. 1. Thus by tiling N a , N b and N c bct unit cells corresponding to fcc or bcc along the a-, b-and c-directions, one can construct a bcc and a fcc supercell which both contain the same number of lattice sites 2N a N b N c . Specifically, the positions of the lattice sites in one of the supercells are given by (n a a, n b b, n c c) and (n a a + a/2, n b b + b/2, n c c + c/2) where n a = 0, 1, . . . , (N a − 1) and similarly for n b and n c .
We have just described how to construct an fcc and a bcc supercell with the same number of lattice sites. For the N V T ensemble we are interested in fcc and bcc phases with the same density ρ, and the supercells should reflect this. It can be shown that a bcc = 2 1/3 ρ −1/3 , a fcc = 2 1/6 ρ −1/3 and c fcc = 2 2/3 ρ −1/3 . Hence for fcc and bcc crystals of equal density it is necessary that a fcc = 2 −1/6 a bcc and c fcc = 2 1/3 a bcc = 2 1/3 c bcc . Therefore, starting from the bcc supercell, if one applies the transformation to the lattice site positions and the dimensions of the supercell, then the end result is an fcc supercell with the same density. This corresponds to stretching the bcc supercell in the cdirection, while simultaneously compressing the supercell in the a-and b-directions to preserve its volume. Taking this idea further, if one makes the aforementioned transformation, but keeps the displacements {u i } of the particles in the supercell unchanged, then one has performed a density-preserving lattice switch from bcc to fcc. Obviously the converse operation is a densitypreserving lattice switch from fcc to bcc. In an N P T ensemble it may be the case that, say, the bcc phase is of a lower density than the fcc phase, in which case a switch which increases the density upon transforming from bcc to fcc, and correspondingly lowers the density upon transforming from fcc to bcc, may yield a more efficient simulation. The above discussion can be easily adapted to treat such non-density-preserving lattice switches.

Results of initial investigations: instability of the bcc phase
Finally we turn to our simulations, the methodology of which closely resembles that of Ref. [7]. Before performing a LSMC simulation to calculate ∆F to a high degree of accuracy, one must optimise the step size used in the particle and volume moves used to generate trial states each time step. To do this it is sufficient to use Metropolis Monte Carlo simulations locked into one of the phases, with a small system size. Furthermore, such simulations act as a 'sanity check' before more accurate simulations are undertaken. It was during such simulations that we noticed that the bcc phase would quickly degenerate into some other -as yet unidentified -metastable phase. The same was not observed to occur for the fcc phase. This is illustrated in Fig. 2, which shows the radial distribution functions (RDFs) at the end of N V T Metropolis Monte Carlo simulations of systems of 250 particles at ρ ≈ 1.1 and β = 3.333 initialised in the bcc and fcc phases, where we are using reduced units as described in Ref. [7]. Each simulation comprised 1 × 10 7 particle moves, and repeats of the simulations yielded very similar RDFs to those shown in the figure. Note that the peaks in the fcc RDF are in excellent agreement with those of the perfect fcc crystal; however, the same cannot be said for the bcc RDF.
The fact that the bcc phase is so short-lived makes it impossible to apply LSMC as it stands to determine the free energy between fcc and bcc in the Lennard-Jones system. A similar problem was described in Ref. [8] for two-dimensional core-softened systems. It should be borne in mind that if one considers any two phases, at least one of them will be metastable, and hence will destabilise given a long enough simulation time; we require that the two phases under consideration do not destabilise before the simulation time required to determine ∆F to the desired accuracy is reached. With regards to the fcc-bcc problem, the following question comes to mind: is there a way in which the system can be kept in the bcc phase for long enough to gather decent statistics relevant to calculating ∆F , or indeed any other property? One might think that a particle move mechanism which constrains the system to remain within the phase under consideration is a valid means of preventing the bcc phase from destabilising. However, 'hard wall' constraints on the particle positions can lead to 'drift' in the centre-of-mass of the simulation, which may invalidate the final results [2]. An alternative approach is to softly 'tether' particles to their lattice sites in the multicanonical simulation through judicious choice of the weight function. A similar idea is used in applying LSMC to fluids [3], and may be worth exploring as a means of addressing instability in solids.