next up previous
Next: About this document ...

A Walk in Phase Space: Solidification into Crystalline and Amorphous States

Joan Adler

Department of Physics, Technion, Haifa, Israel, 32000

The formation of crystalline and amorphous solids is described using a simple glass bead demonstration and a discussion of annealing and rapid quenching on a computer.


This article discusses the use of visualization as a tool for exploring the solidification of a system into crystalline or amorphous states. The demonstrations and the Monte Carlo simulation can be used to help students gain an intuitive understanding of these states without the need for a detailed discussion of statistical mechanics. Alternatively, the material is suitable for a statistical mechanics or a computational physics course after the simulation of Ising models has been discussed.

A crystalline solid exhibits translational invariance so that the atoms arrange themselves in a regular pattern, with a specific lattice spacing between neighboring atoms, and the same number of nearest neighbor atoms for each atom. In the amorphous state the number of neighboring atoms varies, and there is no regular pattern over long distances, although varying amounts of local order may be present. The latter description is similar at first glance to a liquid at a given time, but an amorphous state differs from a liquid in detail. Atoms in the amorphous state do not move far from their equilibrium sites, whereas in a liquid such movement is common. If a liquid is cooled instantaneously (quenched), it will usually go to an amorphous state. In contrast, crystalline states are obtained by slow cooling, with frequent small heatings during the cooling process to remove fault lines. This process is known as annealing.


We start with a square or rectangular clear plastic container. Fishing tackle boxes are perfect. We then obtain 20-100 glass beads of the type used in chemical filters. A size similar to small ball bearings is good. The beads can be obtained very cheaply at chemical stock rooms in most universities. Spread the beads in a single layer that covers about half the base of the box. The box is held at 20-30 degrees to the horizontal and lowered towards the horizontal.

A slow lowering of the box with occasional small shakings mimics the annealing process and gives a perfect crystal (see Fig. 1(a)). Note that the crystalline state in this case is one where each bead has six nearest neighbors and is centered at the sites of a triangular lattice. This lattice is the most efficient way to pack glass beads into a two-dimensional space. Rapid lowering corresponds to a quench and an amorphous structure is obtained (see Fig. 1(b)). Defects such as vacancies and grain boundaries can be induced with a little practice (a vacancy is shown in Fig. 1(c)).

These pictures were made by placing the box on a copy machine. In class, multiple boxes can be used (closed boxes are recommended in this case) or the demonstration can be made on an overhead projector. If using an overhead or copy machine, a coin can be placed to hold the box at a degree or two from horizontal. To use an overhead, clear glass beads are essential; for other uses steel ball bearings would be just as good, but are much more expensive.


We start by formulating a model system. First, we choose the interaction between two atoms. For simplicity, we choose the two-body Lennard-Jones potential V(r) between two atoms a distance r apart to be given by

\begin{displaymath}V(r) = 4\epsilon \bigl [({r \over \sigma})^{12} - ({r \over
\sigma})^6 \bigr] . \eqno(1)\end{displaymath}

The parameter $\epsilon$ sets the depth of the attractive well of the potential and the parameter $\sigma$ sets the range of the hard core repulsion. This potential is accurate for noble gases (except Helium) with the standard example being Argon, and often is used for generic simulations in other systems owing to its relative simplicity. Simulations of real materials, especially metals such as Aluminium, must be made with more realistic potentials than Lennard-Jones.1

Next we choose the number of atoms and their initial positions, and the temperature and the geometry of the space in which they are confined. We also need to think about what type of final state we desire, and the positions and temperature needed in this state. A common initial state is the completely random state. If the desired state is the crystalline state, whose structure may not be known in advance, then we must be very careful that we actually reach equilibrium. The ground state, the lowest energy state at zero temperature, is the triangular lattice if the atoms are allowed to move only in a two-dimensional space. We will consider the canonical ensemble for which the temperature T and number of particles N is fixed, and the probability Ps of a particular microstate s with energy Es is given by2

\begin{displaymath}P_s={1 \over Z}e^{-E_s/kT} . \eqno(2)\end{displaymath}

The function $Z=\sum^M_{s=1} e^{E_s/kT}$ is the partition function, and M is the number of accessible states.

The Monte Carlo method is a way of generating states such that the probability of a state is given by the appropriate distribution for the chosen ensemble. The process of moving from one microstate to another is a stochastic process. This process must sample the states which are most important and obey the principle of detailed balance whereby the ratio of the probability of going from state j to state i to that of moving from i to j is equal to $\exp(-(E_j-E_i)/kT)$.2 The algorithm commonly used to make these moves in the canonical ensemble is the Metropolis algorithm.

The bead systems were constrained by hard walls, but free boundaries with an effective neighbor correction at the edge are used in the simulations. The effective neighbor correction is implemented by using an effective temperature for each atom depending on the atom's current coordination number (number of neighbors within some distance), so that the probability of a change of position is approximately independent of the actual number of neighbors. This scheme facilitates the spatial rearrangement of the particles and so is ideal for modeling the different transformations within a solid. For simulations where spatial reorganization is expected (for example when modeling stresses on lattice mismatched adsorbed systems), this method is superior to the commonly used periodic boundary conditions. However, the method is unsuitable for simulations of a liquid because without walls or an ``infinite'' sample to restrain it, a liquid would decompose into droplets.

A common mistake is to allow the system to become stuck in a local minimum or metastable state and to interpret this local minimum as the ground state. One way to avoid this trap is to monitor a variable such as the mean lattice spacing and to make sure that its mean does not change if the simulation is run twice as long as was originally thought to be sufficient. Another useful precaution, especially for a random starting state, is to ensure that different starting points give the same final result. We also can start at different temperatures to reach the same final result. For long simulations aimed at reaching the ground state, an easy check is to take the final result, make random displacements in each atom's position (unphysical of course), and check that the system returns to the ground state. However, when the amorphous state is the desired final state, we want to end in a metastable state and thus only the route of repeated simulations from different starting points provides a suitable check.

An earlier version of the program described below was given in Ref. 3; the new version has some additional features and is available at our Web site.5 For the present examples we use N=100 atoms. The initial states and temperatures will be varied depending on the desired final state. At each temperature, we find the equilibrium configuration by repeating the basic steps of the Metropolis algorithm until no further changes in some averaged quantity such as the energy are observed. At each step one particle is randomly selected, and a random change is made to its coordinates. The energy before (E0) and after (Ef) the change is compared, and if the energy is decreased, the change is accepted. If the energy is increased, the step is accepted with probability p given by

\begin{displaymath}p = e^{-6(E_f-E_0)/zkT}, \eqno(3)\end{displaymath}

where T is the temperature, z is the coordination number of the particle (number of nearest neighbors), and k is Boltzmann's constant. The factor of six in Eq. (3) is included because the average coordination number in the interior is six. We need to include the local value of z explicitly to make the effective neighbor correction at the edges of the sample. The possibility of accepting a change that gives a local increase in energy enables the system to escape from local minima. Local minima here correspond to fault lines, vacancies, or other defects.

It is helpful to visualize this process while doing the simulations by doing an animated presentation of points in a two-dimensional space. We used PGPLOT,4 which is a simple library of Fortran routines that can be called from C or Fortran and is free for non-commercial use. We have placed our code on the Technion Computational Physics Web site,5 where it may be downloaded for non-commercial use. The code works on UNIX and Linux systems. Also available are C and Java programs.

As in a laboratory experiment, we can grow the desired crystal or amorphous sample by selecting a suitable initial state and temperature, and varying the temperature until we reach the desired final state. To model the solidification process into an amorphous state, we drastically reduce the temperature. An example of such a process is given in Fig. 2, where the initial state was completely random (corresponding to infinite temperature), and the temperature is immediately quenched. In this case we do not want to achieve equilibration at each temperature, because our final state is metastable.

An example of a crystalline state is shown in Fig. 3. All the central atoms are neatly lined up in a triangular structure, but the boundaries are not quite perfect. It is likely that the surface also is not perfectly aligned in a crystal grown in the laboratory. This particular state was begun from a random initial state, but was very slowly cooled to the crystalline ground state. In another attempt to reach the ground state a sample with two differently oriented domains (Fig. 4) was obtained. This system was then reheated and cooled again to obtain the perfect crystal of Fig. 3.

It also is possible to model the processes described above by molecular dynamics. However, molecular dynamics requires more sophisticated programming to solve the differential equations of motion than is required by the present algorithm, and is less directly connected to the principles of statistical mechanics. (The introduction of temperature into molecular dynamics simulations is far more subtle than in the present case, where the introduction of temperature is very natural.) In realistic simulations of solidification it usually is best to do molecular dynamics for some time to model a process with a real dynamics and then to find the local or global minimum with an anneal or a quench using Monte Carlo.

The processes described above are the basis for much current research in computational condensed matter/materials science. Semiconductors, ceramics, Carbon, and polymer systems can be modeled in this way with the proper choice of potentials. Of course, an ab initio quantum mechanical potential that includes the effects of all the electrons is most accurate, but severely limits sample size. When the effects to be modeled are not strongly dependent on the details of the electronic structure, it does not matter which potential is used, if it has parameters which can be adjusted to give the correct lattice spacing and ground state structure. For example, we have seen that Lennard-Jones simulations give results similar to the glass bead system, implying that the hard core repulsion is the most important part of the potential. Carbon must be modeled with a potential capable of describing the insulating diamond structure and conducting graphite and the transformations between them. A comparison between an ab initio and the classical Tersoff potential for Carbon has recently been presented.6 On our Web site5 there are links to Web sites of other groups doing similar calculations. In these calculations, the freedom to move the atoms away from their ground state lattice structure is an essential feature.


1. One situation where crystalline order is not desired as a result of the freezing process is the preparation of ice cream (crystals would spoil the smooth creamy texture). A very fast quench which would give the creamiest taste is sometimes impractical. If we try to freeze ice cream in a freezer, ice crystals start to form, so the ice cream needs to be repeatedly removed and stirred to break up the crystals. This behavior explains why an ice cream machine, which churns the ice cream while it is freezing, gives such desirable ice cream. Suppose that the ice cream at a supermarket has ice crystals in it. Suggest whether it must have melted (perhaps en route to the supermarket), and if so, whether it refroze rapidly or slowly.

2. You could prepare ice cream and justify the time investment and calories with the idea that this preparation is part of your physics homework. You could also try to find other examples of the application of statistical mechanics in food preparation processes. For example, preparing a jello mold is a polymerization phase transition to an amorphous glassy state. Beautiful bubble rafts, sometimes exhibiting crystalline order and illustrating crystal defects can be made while doing the dishes if you are lucky enough to find the right small opening to form regular bubbles from the detergent. Think of other everyday examples.

3. Prepare a glass bead simulation as described above and obtain the different states shown in Fig. 1.

The following problems require downloading the Monte Carlo program5 or writing your own.

4. Determine which annealing schedules lead to crystalline ground states. Some examples are given in Ref. 3. Observe what happens when you reheat from the crystalline phase. Do you see any hysteresis?

5. What happens when you start from the square lattice which is a saddle point structure. Are there differences between annealing and quenching?

6. Think about the differences between crystals and amorphous solids, and how you would determine whether you have a basically crystalline sample with some grain boundaries, or a truly amorphous system. Keep in mind that in nature there is a continuum between small crystallites with grain boundaries and a truly amorphous system with no preferred direction. How difficult is it to make these distinctions which can be seen so easily by direct visualization?

7. Are there observable differences between the structure of a liquid and that of an amorphous solid?


I thank Amihai Silverman for the collaboration leading to the original version of the program described here, and David Saada, Adham Hashibon, and Amit Kanigel for current collaborations in computational condensed matter physics. I acknowledge collaborations with Rafi Kalish and Wayne Kaplan for their encouragement to do simulations to describe their experiments. I acknowledge Alex Zunger for discussions about the importance of moving simulations off-lattice. The support of the German-Israel Foundation during the preparation of this manuscript has been essential.


1.A detailed discussion of interatomic potentials and their applicability to different materials can be found in N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders College, Philadelphia, 1976), p. 398 and H. Gould and J. Tobochnik, An Introduction to Computer Simulation Methods, 2nd ed. (Addison-Wesley, Reading, MA, 1996), p. 214.

2. M. Plischke and B. Beregson, Equilibrium Statistical Physics, 2nd ed. (World Scientific, Singapore, 1994).

3. A. Silverman and J. Adler, ``Animated Simulated Annealing,'' Computers in Physics 6, 277-281 (1992).

4. Information on PGPLOT can be found at$\sim$tjp/pgplot.

5. The Technion Computational Physics home page is at and the page /$\sim$phr76ja/ajp.html contains a summary of links relevant to this article.

6. D. Saada, J. Adler and R. Kalish, ``Computer simulation of damage in diamond due to ion-impact and its annealing'', Phys. Rev. B 59, 6650 (1998).


Figure 1. Structures found from the glass beads demonstration. A perfect crystal (a), an amorphous structure (b), and a crystal with a vacancy (c) is shown. This figure can be viewed in color at$\sim$phr76ja/glassballs.html.

Figure 2. The amorphous state obtained by a rapid quench from a random start to a temperature just above T=0. An animation of this figure can be viewed in color by downloading the file,, and viewing it with Ghostscript.

Figure 3. The crystalline state obtained after a long slow cooling with occasional heatings.

Figure 4. An intermediate state in the cooling process, which after reheating eventually resulted in a picture similar to Fig. 3.

next up previous
Next: About this document ...
Adler Joan