{'fraction_non_alphanumeric': 0.019815913585881636, 'fraction_numerical': 0.0038795070743952533, 'mean_word_length': 4.919180549302117, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nanometric Scale Surface science and the Markov chain Monte Carlo simulation of disordered systems

The Markov chain Monte Carlo method as a statistical mechanics technique for the study of macroscopic systems has furnished the scientific community with great knowledge and advances in the theory of phase transitions. While a number of Monte Carlo models have been proposed for the study of surface growth, these models have not nearly being studied as exhaustively as in the models of magnetic systems, a paradigm of which is the classical model of Ernest Ising. In particular, studies of phase transitions in surface/interface science at nanometric scales are almost non-existent. This article has been written to motivate research in this area of statistical mechanics from the perspective of surface science. In this article we survey the rudiments of the method along with some models of disordered systems such as magnetic systems, material fracture, nano-pattern formation under ion bombardment, and molecular chirality. We performed simulations of these models using the method and obtained results that are in excellent agreement with experimental observations.


INTRODUCTION
The relevance and progress of theories of natural phenomena is strongly rooted in experiments that probe these phenomena under viable and accessible conditions. However, with increasing technological advances the experimental apparatus evolve with increasing sophistication often accompanied by huge operation or acquisition costs that create barriers which scientific enquiries must overcome before experiments designed to investigate them can be approved. This incompleteness in a number of reports, under the general assumption that they are either known or irrelevant, occurs in other areas as well. However, the statistical mechanics theory is more completely explored and developed in the application to magnetic systems and phase transitions [1,2,4,5]. Some of these advances are as a result of what obtains from the application of the theory to seemingly unrelated areas (e.g self-organised criticality in granular mounds, biophysics, etc.), which have also benefited from the reverse feedback of applied results of the fostered advances in statistical theory.
In this article we review the Markov chain Monte Carlo method in general as well as some models of disordered systems to which it has already been applied. The models considered are the Ising model of magnetic systems, a sputter-erosion model for studying pattern formation on material surfaces at nanometre length-scales, a model of material failure, and a model of molecular matching and packing. These models allow for the application of standard techniques of analysis in statistical physics. With recent advances in algorithm designs and combinatorial optimization in theoretical computer science [6,7], the simulation algorithm for these models are being continually developed with increasingly faster and efficient simulations [4] which can help develop other fields and vice-versa in the ever thinning interface between the fields of scientific endeavors.
We start with a review of the basic but fundamental features of a computer simulation of disordered systems based on the theory of statistical physics. We start from standard material and gradually progress to research methods at the level of advanced texts and reference materials. In the context of this established theoretical simulation basis, we then discuss and present results of a few of the existing models

REVIEW OF THE MARKOV CHAIN MONTE CARLO METHOD
The macroscopic or thermodynamic properties of a system can be calculated as a function of its microscopic properties through the relation (see e.g. Ref. [5,9] Depending on the nature of the interaction of the system with its environment there may be one or more of energy exchange, matter exchange, or exchange of spatial extent, which determines the statistical ensemble, i.e. assembly of microstates subject to one or more of these constraints, within which the system evolves and is simulated.
In any such scenario the equilibrium state of the system is any of a large number of equally probable microstates, whereas the probability of a microstate varies between equilibrium states. In fact the number of microstates is infinite if the phase space points form a continuum.

MARKOV CHAIN MC SIMULATION OF DISORDERED SYSTEMS
For reversible processes carried out on a canonical ensemble where and the number of particles are constant, and the number of microstates with energy at temperature is given by The relationship between and is that between the average of a data set and the elements of the data set; as we shall see below the system must be allowed to equilibrate at a particular temperature before measurements are made, nevertheless, the energy fluctuates about an average value after the equilibration time such that ; where is the number of microstates with (possibly different) energies through which the system evolved in the measurement process. In the canonical ensemble is given by and is given by (7). is not necessarily normalized to have unity integral or sum over the sampled region but is proportional to a probability. Now is a conditional probability, i.e. the probability of accessing the microstate given that the system is in the microstate So are transition probabilities of transition from microstate to . For us to obtain the partition function we would have to first sample all microstates and evaluate the sum of the distribution functions which is, however, unnecessary since we must have that is, The condition (11) is known as detailed balance. It allows us to calculate transition probabilities in terms of relative transition probabilities between successive microstates without actually calculating the partition function. For instance, in the canonical ensemble (7) and (10) imply where is the energy difference in a transition from the microstate with energy to the next microstate with energy in the Markov chain, and is the relative conditional probability that the system will evolve to microstate from the microstate .
According to (12) the system will certainly change into the microstate from if since in this case , whereas it will likely remain in the microstate if since in this case and as . This ensures that the vast phase space is "importance" sampled in such a way that the system evolves in the few most probable microstates in the phase space in the course of the simulation and reduces the rigor and inefficiency of having to go over or uniformly sample ("simple" sampling) the entire phase space in the calculation of average quantities such as, e.g.
In the calculation of (13) most of the Boltzmann weights are vanishingly small, corresponding to the very large number of microstates in which the system is unlikely to be found, and therefore contribute very little to the average value.
Note that, depending on the application, and may have entirely different physical meaning.

Simulation Algorithm
Based on this theory, a Markov chain MC simulation of the evolution of the microstates of a system in statistical mechanics can be performed by using the following algorithm, known as the Metropolis algorithm: Generate a trial configuration for the initial microstate of the -particle system and obtain the next microstate in the Markov chain by tweaking the coordinate of a single (  This algorithm therefore ensures that the system always tend towards microstates of lower energy than the previous one, but unlike the so-called "greedy algorithms" [10] that never accept a higher energy option, it also allows the system to sometimes (though with very low probability) accept a move to a higher energy microstate, wherein lies its power and vast applications to diverse optimization (minimization/maximization) problems. For instance, if higher energy moves are always rejected then there will be no possibility for the system to move out of any of the local minima in Fig. 1 and the system would therefore not equilibrate to the global minimum. However, as the temperature decreases the probability of the system moving out of a local minimum, if trapped there at low , also decrease and at such low it may be difficult to attain a microstate of global energy minimum; depending on the energy landscape of the system.
It is best to run independent simulations with different initial conditions, e.g.
one with a cold start ( ) and another with a hot start ( ), so as to avoid inaccurate results due to local minima pitfalls. For the cold start a starting configuration or microstate that is known to have the lowest energy is used, while for the hot start a configuration that corresponds to the highest energy is used. The result of the independent simulations must be the same irrespective of the starting configuration, otherwise one or both must have ended in local minimum and further checks are necessary.
Since the system equilibrates to a "global" energy minimum in which it then spends the rest of its time, the probability is time dependent , i.e. and its time dependence is governed by the master equation where is the probability of the system being in the microstate at time , and is the rate of transition from microstate to . In equilibrium and .
Obviously the efficiency of this algorithm, at the so-called "importance sampling" of the phase space region in which the system spends most of its time, will depend on the energy landscape of the system; the algorithm will be more efficient and faster for a system with the energy landscape of Fig. 1(a) than that of Fig. 1(b). As a result, a number of cluster algorithms have been proposed to increase the simulation speed and efficiency at "importance sampling" by avoiding the slow changes of state that occur with time in certain stages of the simulation (e.g. critical slowing down, frustration). Such algorithms consider microstate transitions involving clusters of particles instead of single particles [11 -14].
Since the system may end up in local minima it is better to let the simulation run for a while before taking measurements so as to allow the system to equilibrate; the time required for the system to equilibrate is known as the equilibration time, which depends on the energy landscape of the system. The configurations, or microstates, are usually stored as arrays whose elements are the single-particle states corresponding to each particle of the system.
The energy of a microstate depends on the interactions between the particles and is the sum of the energies of these interactions. The particles, being indistinguishable (e.g. fermions or bosons), must have identical surroundings and interactions in such a way that the boundary particles are not distinct from the inner ones. This is achieved by imposing periodic boundary conditions on the system, thus transforming the -dimensional system into a -dimensional torus, so that the particles at one end of the system sees the particles at the other end as nearest neighbors. The storage requirements for specific problems or the geometrical nature of certain systems may require other boundary conditions like screw-periodic, antiperiodic, mean-field, or free-edge etc., boundary conditions. Other lattice structures may also be more appropriate for specific problems [1,2].

Uniform Deviates
The recipe for simulation is incomplete without (pseudo) random numbers [10]. We can not choose in a deterministic version because that would be unrealistic, it would either imply a prior knowledge of the phase trajectory followed by nature or that we are the ones dictating what course nature would take which is impossible for physics in the subatomic realm. Hence, is chosen randomly by using floating-point uniform deviates such that each has the same likelihood of being picked as any other. Our knowledge of the probability of a trial microstate is also not enough because we need to know when to accept or reject it according to its probability. Again, a floating-point uniform deviate is used.
For instance, if the probability of an event is 0.1 then the event will occur once out of every 10 attempts, as the number of attempts tends to infinity. To simulate this uniform deviates are used. The random numbers are assumed to be uniformly spread over an interval (0 to 1 for floating-point random numbers) such that each random number has an equal probability of being chosen as any other one.

Finite Size Effects
The number of particles (system sizes) that may be considered in a simulation is much less than the number of particles in a typical macroscopic system, i.e. the simulated system size whereas in the thermodynamic limit we should have . Typically, real particles whereas simulation . For instance, a system with only two possible single-particle states will have microstates or configurations while the simulation phase space comprises of microstates.
However, use of small still gives very accurate results when compared with real systems for which provided the simulation is for temperatures far from the critical temperature at which a phase transition occurs. For temperatures close to and above the critical temperature the effects of using small , known as finite size effects, become very strong and the simulation results become inaccurate; the inaccuracy increases with decreasing . Hence, for studies involving phase transitions, or for simulations at temperatures above , it is important to do a finite size scaling of the results to obtain the actual results of a real thermodynamic system.
The form of the finite size scaling will depend on the set of thermodynamic quantities being measured and their scaling relations (1) (2) (3). It also depends on the order parameter, which is a quantity that has a value of zero on one side of the phase transition and a nonzero value on the other side of the phase transition; i.e. which is like a step function whose step occurs at the phase transition.
The origin of the finite size effects is in the divergence of the correlation length (the length scale of the fluctuations in the order parameter) at the critical temperature as a result of which it is impossible to get accurate results of the thermodynamic properties of the system at from any simulation performed on a finite lattice. Since the lattice size is less than the value of at all the lattice points are correlated and one is unable to obtain the scenario of linearly independent singleparticle states.

SIMULATION RESULTS.
In this section, we present some general simulation models which are typical examples of the successes of the Markov chain Monte Carlo approach to the simulation of the complex, disordered, systems frequently encountered in nature which are of great technological importance.

Magnetism and the Ising Model
Magnetism is one of the most studied phenomena in materials. It has been a study of active and intense research interest for centuries [15,16]. Any material is where is the unit of the spin-spin exchange interaction and is often chosen to be , in which case the lowest energy state will tend to have neighbouring spins aligned and the low temperature ground state will be a ferromagnet (all spins aligned), or in which case the lowest energy state will tend to have neighboring spins opposed in orientation and the low temperature ground state will be an anti-ferromagnet (alternating spins). denotes a nearest neighbor pair, and is the unit of the spinfield coupling.
The spin operators (actually for each spin, assuming the field is oriented along the -axis so that ) have eigenvalues where is the Dirac action constant. In simulations the eigenvalues of the spin operator are chosen to be . Thus the number of single-particle microstates of the system is reduced to only two; for up and down orientations of the intrinsic spin.
When the exchange interaction energy varies from bond to bond then we have a spin-glass system which is widely studied with a variety of interesting ramifications [17].

Model.
A more recent observed phenomenon in surface science, than the magnetism of condensed matter described above, is the formation of self-organized nano-patterns on material surfaces when bombarded by a beam of energetic particles, with energies in the range of keV [8, 18 -23]. The surface morphology arises as a result of the interplay between the roughening processes of stochastic removal of surface particles by ion-sputtering, and the smoothening processes of diffusion of surface particles or radiation-induced viscous flow. The roughening process creates instability of the surface against further perturbations by the eroding ions such that the sputter yield becomes curvature dependent with the result that surface depressions are eroded in preference to surface protrusions. According to the continuum theory, it is this curvature dependence that leads to the formation of ripples of wavelength , where is the surface diffusivity, and is a surface tension coefficient. The wavelength of the ripples is on the nanometer scale and they are oriented along the direction with the largest [21].
The experimentally observed ripple orientation is parallel to the projection of the ion beam direction onto the surface plane for large incidence angles (close to grazing incidence), and perpendicular to the projection for small incidence angles; except for metallic surface with anisotropic diffusion in which case the ripple orientation is perpendicular to a crystallographic direction (i.e. the one favored for diffusion) at small incidence angle.
Hartmann, Kree, Geyer, and Koelbel (HKGK) proposed a discrete solid-onsolid Monte Carlo model [24] for the simulation of surface modification and evolution by ion-bombardment which, in contrast to the Ising model above, simulates the surface morphology indirectly by resorting to the effect of collision cascades on the surface evolution without the direct need for a surface Hamiltonian. These collision cascades are set off in the near sub-surface layer by the impinging ion and are responsible for the redistribution of the energy of the impinging ion to the surface particles. The energy received by a surface particle at position due to the arrival of the ion somewhere in the vicinity of this surface particle is assumed to be of the Gaussian form [25]: where and are the widths of the collision cascade ellipsoid parallel and perpendicular to the ion beam direction, respectively. A sputtering process is simulated by eroding surface particles in the vicinity of the impact position of the ion with probabilities proportional to [26 -28], and a diffusion or smoothing process is simulated by using any of the available diffusion models (e.g. [29,30]), depending on the experimental hopping rates for the material concerned.
We performed extensive simulations using this model and the surface profiles obtained from a typical simulation are shown in Fig. 2

Molecular Chirality and the YNL Matching Algorithm: Simulation Results.
The problem of measuring the extent to which an object fails to coincide with or match its mirror image [31 -39] when optimally superimposed on the mirror image was recently tackled in a direct way through simulations of the matching process using the Markov chain MC method [40] (YNL algorithm). The measure of chirality or degree of matching of the objects (an object and its image or two different objects as in molecular similarity or facial resemblance) that have been optimally superimposed by the YNL algorithm is based on the concept of the Hausdorff distance between sets as proposed by Buda and Mislow [33].
This direct approach entails the simulation of the matching process which includes a large number of rotations and translations such that the Hausdorff distance between the objects being matched is as small as possible. if an object is achiral or if two objects match, and if an object is chiral or two objects do not match. The challenge is that one may get very close to the global minimum of only to be pushed into a completely different region of the landscape with the slightest of translations in the matching process. But the YNL Monte Carlo matching algorithm [40] always arrive at the optimal match. Hence, this algorithm can be of immense benefit to the handling of harmful enantiomorphs of useful chiral molecules in drug design, to database searches in the area of sequence alignment and queries on other databases, and to the understanding of the mechanism of chirality transfer from a chiral dopant to an achiral host.
The YNL algorithm searches for the optimal superimposition of an object and its mirror image , or second object, by performing a series of random rotations and translations of keeping stationary. The simulation is started with a random orientation and translation of which represents the initial state or configuration. The matching process is iterated down over the 'temperature', thus annealing the process and ensuring that as the system cools, the rotation and translation reduces in magnitude so that at high 'temperature' the system moves freely from one region of the landscape to another in the search for the region that contains the global minimum, and much less freely when this region is supposed to have been found (at low 'temperature'). During each iteration stage over the temperature, it is kept fixed while Monte Carlo steps of the matching process are performed. In each MC step the state of the system, which is the orientation and position of relative to , is changed by rotating and translating and the move is accepted according to the standard MC transition acceptance rules, tailored towards acquiring minimum or . provides a measure of the chirality of the object through the chirality measure of Buda et al. [33] where is the minimum Hausdorff distance between and over all positions and orientations of and , is the diameter of . The results of our simulations as applied to the multidimensional search for the most chiral tetrahedral plane shapes with D 2 , C 2 , and C 1 symmetries are presented in Table 1. As can be seen from the table, our results are in excellent agreement with the results of Buda, Auf der Heyde, and Mislow [32,33] who obtained their results by using another, albeit more tedious, numerical optimization method.

Material Fracture and the Fiber Bundle Model
Fracture in heterogeneous materials is a complex physical problem which has been a focus of research interest for a long time [41 -43]. Fiber bundle models (FBMs) form a fundamental class of approaches to the fracture problem through their capture of the essential features of material breakdown. They are models of materials as composites of fibers such that the breakdown of the material is as a result of the propagation of fracture among the failing fibers. One such model which provides a deep understanding of the intrinsic nature of the fracture process is the dynamic fiber bundle model [44]. According to the dynamic FBM the material breaks down as a result of the fatigue of its fibers over time. That is, when the material is stressed or loaded the fibers share the load and get progressively stressed with time such that even if the external load or source of stress is removed the acquired fatigue in the fibers remain and are incremented with the introduction of a new source of stress until the individual fibers reach their threshold and fail, leading to the overall failure of the material with time. This differs from the static FBM [45] which attributes the failure of the individual fibers to their quasi-static loading due to the gradual increase in the external load.
The dynamic FBM includes a normalised combined load sharing rule [44] is the distance between an active fiber and a failed fiber , denotes the set of active fibers, and is a variable parameter that controls the effective range of interaction among the fibers. This combines the local load sharing rule ( ) in which the load borne by a failed fiber is shared equally among the four nearest neighbour fibers, the global load sharing rule ( ) in which the load of a failed fiber is shared among all existing fibers, and the variable range load sharing rule ( ) in which the load of a failed fiber is shared beyond the nearest neighbours but not globally, with the range depending on .
The probability of the breaking of a fiber in one sweep of the lattice in the time interval is given by is the Weibull index ( ) which gives the degree of heterogeneity of the system; as increases the material becomes more homogeneous.
is the stress on the fiber at time . The results of our simulations are presented in Fig. 3. Our results indicate that the failure time of the material decreases with increasing range of interaction between the composites of the material and this decrease is precipitous when the interaction range is just above next nearest neighbour interaction. Also, our results indicate that the time to failure is delayed if the material is made less heterogeneous; that is, of fewer composites. Again, these results are in excellent agreement with experimental observations [46,47]

Conclusions
This article reviews the modern theory and implementation of computer simulations disordered condensed matter systems and fluids, based on the statistical mechanics approximation method of Markov chain Monte Carlo techniques. We started with the fundamentals without recourse to any system or problem for a generalized treatment and then discussed a number of recent applications including Ising model or spin-glass model of magnetic systems, a solid-on-solid model of nanostructure formation on surfaces driven by ion-bombardment, statistical MC models of material fracture, and a MC model of chirality which is applicable as a matching algorithm. Basic and typical results of our simulations of the different complex systems are presented.
It is hoped that the in-depth treatment of the rudiments of the method, and the excellent agreement of its results in diverse applications would motivate and stimulate its application to the solution of still unsolved statistical physics problems from the perspective of the richer energy landscapes to be encountered in nanometric surface science.