The Markov Chain Monte Carlo method was first conceived by Stanislaw Ulam, a physicist working on mathematical physics problems at Los Alamos as part of the Manhattan Project, he was in the hospital playing solitaire when he decided to try and calculate the probability of dealing a winnable game of solitaire[2,4]. After the traditional combinatorial methods had failed him, he’d decided that perhaps it was easier to deal multiple games and observe the number of wins in order to determine the probability of dealing a winnable game[2]. By dealing enough games of solitaire and observing the outcome out each game, it would be possible to construct a distribution of outcomes of the random samples from which the probability of dealing a winnable game can be estimated. Ulam then began to wonder how this method of random sampling could be used to solve various problems that he and his team were faced with at Los Alamos.
Although random sampling was a widely known statistical tool, Ulam’s ideas were revolutionary in that he wanted to utilize computers to simulate the process of random sampling. Ulam began thinking about how the continuous process of change for a system could be simulated in a step by step process with random sampling[2]. Jon Von Neumann heard about Ulam’s ideas and thought to implement it and obtain multiplication rates of neutrons through fission processes in the various designs of the bomb. In particular, they wanted to investigate how long a neutron would travel before it collided with another particle and multiplied through the fission process[2,4]. Since their method left certain sampling outcomes up to chance, they gave the method a name that reflected this characteristic, the Monte Carlo method. The method proved to be successful and provided the scientists with reliable estimates of neutron multiplication rates and allowed the scientists to monitor the path taken by a neutron as it travelled[2]. Sometime afterwards, Ulam and Hastings developed the more widely recognized Markov Chain Monte Carlo that produces random samples from a Distribution of interest by allowing the system to sample all possible outcomes by utilizing a Markov chain and adding in extra accept/reject step to ensure that all possible states, no matter how unlikely can be sampled. In order to use such a powerful method, it is necessary to first understand the Monte Carlo method and Markov Chains separately before taking on the Marko Chain Monte Carlo.
Monte Carlo is a method of obtaining an estimate of properties of a Probability Distribution by obtaining random samples directly from the distribution of interest[3]. The set of random samples taken will yield their own distribution that acts as an estimate of the target Probability Distribution; this estimate can be used to study the behavior of the target distribution. The greater the number of samples taken from target distribution, the better the estimates that can be made by the sample distribution. An example of how one might employ Monte Carlo methods would be flipping a fair coin multiple times in order to obtain an estimate of the Probability distribution of outcomes for the coin. By recording the outcome of each flip, an estimate of the distribution of outcomes for the coin can be observed; as the number of flips is increased, the better the estimate of the true probability distribution will be.
If the Probability distribution function (PDF) is known, then it is possible to simulate a Monte Carlo and produce random samples from the Probability distribution without performing any experiments or obtaining real samples from the distribution. The PDF is a function that describes the probability of observing any particular outcome from a set of possible outcomes; if it is known, it is possible to assign probabilities to each outcome and perform a simulation that will generate samples from the target PDF.
Figure 1
In the example above, a Monte Carlo simulation was performed to produce random samples from the Standard Normal Distribution. As can be seen in the Figure 1 above, the sample distribution becomes increasingly similar to the normal distribution as the number of steps increases. In the Monte Carlo simulation that was performed by Ulam and Von Neumann, the region of values between zero and one was divided into distinct sections that resembled the probability of outcomes[2]. At each decision step, a random number was randomly generated from a uniform distribution and its value would determine what event occurred[2]. By using a Markov Chain, the history of events illuminated the path taken by a single neutron.
A Markov Chain is a series of transitions that a system undergoes in which the present state of the system is a result of the previous state of the system. The system begins in an initial state and can start in any initial state. The system will then transition to a new state, chosen at random, for a prescribed number of steps, the random choice can be completely described by a transition matrix; a matrix that contains the complete set of probabilities of the system transitioning into any state given its present.
Figure 2
The system in Figure 2 above can only assume two states, 0 and 1. The initial state is sampled at random from the initial state distribution shown in Figure 2B. Once the initial state is sampled, all following states are sampled using the transition matrix shown in Figure 2C. The next state of the Markov Chain is chosen as follows: if the present state of the zero is 0, then the probabilities in row 0 are used to sample the next state of the system. This way, the chances of the system assuming a particular state in the next step depends on the present state of the system.
An interesting trait of Markov chains is that if they are allowed to run for a large number of steps, the distribution of observed states will achieve a stationary distribution. The stationary distribution describes how many times each of the possible states was present over the history of the Markov Chain.
Figure 3
For example, if a state that can assume three possible states (A, B, and C), the probability distribution of assumed states for the path taken by the Markov chain will vary for short paths. However, if the number of steps is increased and a longer path is taken, the distribution of observed states will eventually stop changing and assume a stationary state as is seen in the Figure 3. The stationary distribution can be determined from the transition matrix of the Markov Chain using methods of linear algebra, however, these methods are beyond the scope of this discussion.
The Markov Chain Monte Carlo is a Markov chain whose stationary distribution matches that of a target distribution and whose history provides random samples from that target distribution. The MCMC algorithm uses an accept/reject step to accomplish the sampling. However, the most useful aspect of MCMC is that only the unnormalized PDF of the target distribution is necessary to perform MCMC. This makes MCMC very useful in drawing samples from distributions that or difficult or impossible to sample from directly[3].
The Markov chain of MCMC allows the system to assume various states. However, rather than simply sampling a new state from a transition matrix, a new state will be proposed and will be accepted or rejected based on the outcome of a random event. The outcome of each step represents a sample taken from the desired distribution, the accept/reject step ensures that every outcome of the sample space is sampled. Since the stationary distribution of the MCMC is the target distribution, an increase in the number steps taken by the Markov chain of MCMC will yield a more reliable estimate of that distribution.
The data that is obtained from an MCMC method allows for one to study and observe behaviors and important characteristics of the target distribution. There also exist a wide variety of MCMC methods in which the process for proposing transitions vary yielding unique types of paths for the Markov Chain. For this reason, MCMC methods enjoy a wide range of application within all academic fields of study. In this post, we’ll be employing MCMC to construct our own physical-numerical model that will model the formation of a liquid water bridge between two surfaces.
In this investigation, a physical-numerical model will be developed, based on the Metropolis-Hastings algorithm, for the formation of a water bridge between two surfaces that are very close together in an effort to understand how pressure changes and influences the formation of the water bridge in a nanoconfined space. Just as Ulam and Von Neumann did in the 40’s, we’ll be simulating the process of change of a system over time by using our knowledge of the equations of state and energy calculations to produce a series steps that, when viewed, will show the eventual formation of a water bridge. The information that can be derived from such a model has been used1 in the design of Atom Force Microscope (AFM) tips utilized in Dip-Pen nanolithography(DPN). DPN is performed in an environment that is open to the atmosphere, meaning, there are water molecules in the gaseous phase that are capable of forming interactions with the microscopic tip. In DPN, the AFM tip and surface are both very close together which allows for condensation to occur at a much lower pressure than the bulk saturation value1. From the condensation that forms on each surface, a water bridge can easily form between the two which gives rise to capillary forces which influences the resolution of images produced in DPN[1].
In order to develop an appropriate model that provides an accurate description of a system, or at least a decent working model, one must first think about the system being simulated and how certain changes affect the system over time. The system that our model will simulate in this experiment is that of water molecules in a gaseous phase transitioning into the liquid phase through the application of pressure by the reduction of space between the two surfaces. From this description of the system, and from phase diagrams of the system, the condensation of water seems to be a continuous and simple process, however, what occurs on the molecular level is more chaotic.
When the water molecules reach the space between two surfaces, the pressure increases between the two surfaces since the volume of space is decreased. The reduction of space between the two surfaces increases the likelihood of collisions between the water molecules with the surface as well as other water molecules. With the reduction of space and increase in collisions between water molecules, the time the molecules spend near each other, or the surface, increases resulting which allows for the formation of strong intermolecular interactions. As the molecules begin to condense, the energy of the system achieves a lower energy state since the formation of the Hydrogen bond interactions release energy into the surroundings. The formation of longer-lived interactions is a sign of condensation, once the process begins more water molecules will aggregate until a bridge of water molecules between the two surfaces is formed.
Although the molecular description of the formation of the water bridge seems complex and chaotic, it can be simplified greatly by imagining the whole process as a series of still shots compiled together in the form of a movie. Each still shot shows the same system in a slightly varied configuration. Some of the water molecules have moved to slightly different positions while others may have stayed put.
Now that we’ve picked a system to break down, we’ll need to identify a way to simulate the changing states of the system in a way that obeys behavior predicted by the equation of state that defines how the system behaves at equilibrium. To simulate the formation of a water bridge, we must think about the interactions that occur between molecules in addition to how their interactions result in changes of the overall state of the system. Most importantly, we must find a way to translate the continuous system changes described by the equation of state into a system whose changes can be tracked through a Markov Chain and how the accept/reject step can be integrated into the process of transition of the system from one state to another.
For the development of this model a lattice will be used to represent our system since it is easy to construct and manipulate in coding space. A lattice is any space divided up into to smaller spaces of a particular size that can be occupied by a set number of items making them useful tool in keeping track of positions of molecules within the space constraints of the model. In our model the lattice will be represented using a 2-d matrix in which the lattice spaces will represent a molecular diameter. Each position within the matrix site will be capable of holding only a single molecule and will thus be described as being occupied or empty. The occupancy of each position in the 2-dimensional will be described using numerical values since they can be more readily recognized and acted upon when running the algorithm; 0 = empty site, 1 = occupied. In reality, the system will change states by the random motion and collisions of the water molecules as they move through space and collide with each other, the DPN tip, and solid surface. As the water particles collide with the DPN tip and the solid surface, some will condense on these two surfaces; over time, the volume of condensate between the two surfaces increases and eventually gives rise to a bridge of water molecules that connects the DPN tip and the solid surface. Our lattice will emulate the random motion of particles by the randomly changing positions of the occupied sites and unoccupied sites within the lattice.
In order for the random swaps to move towards the formation of the water bridge in the lattice model, the energy of interaction between the molecules and the driving force of the change in energy of the system to a lower energy state must be accounted for. The Ising model will be used to calculate the energy of interaction of the particles and the change in energy of the system. The lattice model is incapable of simulating collisions between particles, however, it is possible to simulate the aftermath of a collision in the form of intermolecular interactions as two particles sitting next to each other, “nearest neighbor” pairs. Using the “nearest neighbor” interactions, the energy of the current state of the system can be determined allowing for energy calculations to be made and taken into account allowing for a model that accurately moves towards a lower energy state.
Figure 4
With a computational representation of our image decided and a loose outline for how the MCMC will be performed, a general flow chart can be written for the algorithm to organize what types of inputs will be necessary for each function and their respective outputs. The flowchart shown in Figure 3 shows the general outline steps that must take place in order to produce an MCMC that will sample the system, in particular, the physical chemical model built in this investigation will follow that of the Metropolis-Hastings sampler (MH). The system will be spend more time in states that are more probable and less time in states that are not so likely as described by the shape of the PDF. In the MH algorithm, if the proposed state has a higher probability in the posterior distribution then the proposed state should automatically be accepted, else, the decision to accept or reject the proposed state will pass to a specially designed accept/reject step3. In the case of physical-chemical model that will be built for the formation of the water bridge, the energy of each state will be used to determine its likelihood of occurrence which will ultimately bias the accept/reject step of the model to accepting lower energy states.
Based on the flow chart shown previously, there are various individual steps that can be written as unique functions that will take in prescribed inputs and generate the desired output. However, within each of these functions there will likely be smaller nested functions that are just as important as the main functions outlined in the flow chart. Once these functions are written and put together in the larger scheme of the MH algorithm, a fully functional model emerges as is shown in diagram one of the appendix. However, before the algorithm can be fully achieved, all necessary functions must be written before they can be pieced together and the MH algorithm achieved. Based on the flowchart for the algorithm, there are at least four obvious functions that must be written: a lattice generator, energy calculator, a swap function that will be capable of deciding whether to accept or reject the proposed state.
The first function in of the algorithm must generate the a randomly configured 2-dimensional lattice with N number of particles, lattice dimensions of (x,y), surfaces for the particles to adhere to, as well as updatable lists of positions of particles and empty space after each step. Since python has special modules that can be used to make the generation of a lattice simple, the only inputs required for this function will be the desired number of particles and the dimensions of the lattice.
Figure 5
Since, our model does not account for the differences in identity of the surfaces and water molecules, some creativity was involved in ensuring that the appropriate number of water molecules were present in addition to the surfaces that line the top and bottom of the lattice; the number of water molecules and surface particles were combined into the input value N such that the desired output lattice would be generated as shown in the output column.
Once the lattice is generated, if follows that the next function to be written is that which calculates the energy of system in its present state. Since the Ising model is being used to calculate the energy of the system, the energy calculation will have within it a smaller function which counts the “nearest neighbor” pairs. Once the nearest neighbors are calculated, the energy of the system can then be calculated as outlined previously.
Figure 6
The function count_nbrs peforms the task of counting all nearest neighbors of each water molecule listed in the “Particle index” list within the input lattice. The function locates the position of the particle and looks above, below, left, and right of the particle to determine if there are any particles present and returns the number of “nearest neighbors” around the particle. Inside of the energycalc function, the nearest neighbors are counted and summed for over all positions in the particle index. This summed value of the nearest numbers is multiplied by 0.5 in order to correct for double counting of pairs. The neighbors are then multiplied by beta, whose negative value places the energy of the Ising model in the context of thermodynamics in which favorable energy states have negative values.
The last function to be written will be called the swap function, a function capable of proposing new states and deciding whether to accept or reject the proposed state. This particular function works by a water molecule at random from the list of occupied sites and selects an unoccupied site at random from the lattice. Since the swaps are proposed by moving a single molecule at a time, the difference in energy between the current and proposed states is the difference in “nearest neighbor” interactions at the current position of the molecule and of the proposed new location. To reduce the amount computational resources consumed, the runtime of the simulation, the difference in energy between the current and proposed states will be calculated by the difference in nearest neighbor interactions at the current and proposed locations.
Figure 7
Once the energy difference between the states is known, the new state will either be accepted or rejected based on the difference in energy when compared to a random number. If the new state is accepted, the list of occupied positions will be updated with the particle’s new location and the lattice will be updated to reflect the new state.
With all necessary pieces of the MH-algorithm written out, they can then be combined into one single function that will perform the MCMC for our system. Since this function has nested within it all of the previously discussed functions, all the inputs required for those functions must be defined so that the algorithm can run to completion. In addition to simply taking in the inputs, generating a lattice, and running the swap function multiple times, the function shown below produces a very helpful output, a histogram which contains the probability of a particular location within the lattice being occupied.
Figure 8
This histogram, coupled with the trajectory of the Markov Chain, can provide insight into the formation of a water bridge when the surfaces are close together and the volume between the two surfaces is held constant. In order to study how pressure influences the formation of a water bridge as the pressure of the system changes, the MH algorithm must be performed multiple times at various lengths such that the distance between the surfaces varies.
In the simulation, the number of steps that the simulation allowed to run is akin to the amount time that is allowed to pass, therefore, the greater the number of steps will likely increase the likelihood of us seeing the formation of the water bridge. If the first frame to the last frame of the MCMC simulation are compared, it can be seen that the particles will have moved much closer together indicating that the particles have congregated, thus, indicating a transition a transition to the liquid phase. Figure 9 below shows the MH algorithm run multiple times with varying rows of space between the two surfaces, but no increase in the number of particles. Since the number of particles was not increased from one run to another, the more spacious systems do not show the formation of a water bridge but do show aggregation near the surfaces which is a promising sign of possible condensation.
Figure 9
Currently, the model does not take temperature into account but by including the temperature into the system the value of Beta would be scaled by the temperature and would impact the probability of the formation of a water bridge. Higher temperatures would make it more difficult to form a water bridge since the particles move with greater velocities, they would be less likely to stay close to each other after collisions take place. At lower temperatures, the probability of observing a water bridge would be higher since the particles move very slowly and are likely to stay next to each other after a collision. However, the variation in bridge formation observed depends on the dimensions of the lattice and the number of particles inside the lattice.
Based on the preliminary results of this model, it does provide qualitative data that meets the initial expectations for the model. The model does show the condensation of water molecules near the solid surface, however, the observation of a water bridge is not very clear and would require a more thorough statistical investigation. Additionally, the model still lacks many features of the system it is trying to simulate such as long range interactions between water molecules. Although this model is quite simple, it provides a starting point towards the building of more complex models as well as useful thermodynamic data can be harnessed an analyzed if a more mathematical treatment is applied.
- Jang, Joonkyung, Schatz, George C. “Lattice Gas Monte Carlo Simulation of Capillary Forces in Atomic Force Microscopy.” Journal of Adhesion Science and Technology. 2010. 24 P 2429-2451.
- Eckhardt, Rogers. “Stan Ulam, Jon Von Neumann, and the monte carlo method.” Los Alamos Science Special issue 1987.
- Ravenzwaaiji, Don Van, Cassey, Pete, and Scott D. Brown. “A simple introduction to Markov Chain Monte Carlo sampling.” Psychon. Bull. Rev. 2018. 25 p 143-154. DOI: 10.3758/s13423-016-1015-8.
- Robert, Christian, and George Casella. “A short history of markov chain monte carlo: subjective recollections from incomplete data.” Statistical Science. 2011. 26.1 p102-115.