The fast algorithm for the numerical simulation of crystal growth.

by Dr. Leonid Sakharov

There are principle methods of algorithm for the numerical simulation of crystal growth these permit dramatically increase its performance by means of pre-calculating all intermediate numbers and presenting them as integers.

The fast algorithm for numerical simulation of crystal growth is the code implementation of the the model of surface of growing crystal described in a page " Physical background of the model for numerical simulation of crystal growth". The model is based on the premise that surface of the growing crystal can be presented as an two dimensional integer matrix L and the scientific background for calculating probabilities of each molecule (that is represented by one element of the matrix L[i,j]) to change positions at the surface of the crystal on the base formulas derived from thermodynamic characteristics of the phase transition and crystallographic parameters of the crystal.

A software implementation of the model for the numerical simulation of the crystal growth have to be as much efficient as possible to permit collection large enough set of numerical experiments for multivariate analysis of results. In our version implemented in LeoMonteCrystal that is written in Visual C++ environment the major efforts was concentrated on the minimization part of operations with float type of variables (these have decimal points in them).

In spite that formulas for probabilities for molecule movement unavoidable demands float point operations like exponent function there is a way to avoid letting them inside loops of simulation itself by following means:

• - parts of surface associated with different directions of cuboid are represents with integer values by scaling whole molecule surface to big enough integer number Sm (to minimize error of rounding) and rounding values of the parts to integers;
• - represent probabilities for molecule to stay in place or to move ahead by integer numbers scaled between 0 and Rmax that is maximum value that generator of random numbers cam produce.
• - calculate two integer two-dimensional matrices with size Sm×Sm there probabilities for molecule to stay in place or move ahead are prevented by integer numbers in the interval 0 and Rmax for all possible combinations of the molecule surrounding on the surface. Note that size of matrices has to be Sm×Sm not just Sm because by our model of quasi-independent behavior of molecule one need to take into account not only current surrounding of the molecule but also its virtual surrounding if the molecule step one move ahead.
• - use integer numbers characterizing molecule surrounding as an indexes in arrays described above for updating probabilities of each molecule to stay or move ahead.

The scheme of fast algorithm for calculation of probabilities of one a molecule movements.

One of most common operation during numerical simulation of crystal growth is defining probabilities for a molecule at given position surface to change its value that it is equivalent movement of crystal surface at this position. Each element of matrix L[i,j] have probabilities to add 1, lose 1 or retain the same value. These probabilities are defined for each time step with two integer matrices with size identical to matrix L: matrix Po[Nb, Nc] contains probabilities for molecule to stay at the same place and matrix P1[Nb, Nc] contains probabilities to move ahead. Integer values of elements Po and P1 represent probability normalized in range of values from 0 to the  maximum integer number Rmax that is specific for given generator of random numbers used in software. For implementation of the algorithm in LeoMonteCrystal this value is Rmax = 0xffffffff = 4294967295.

The core procedure of simulation is repeated in infinite loop can be described as following:

During update of matrices of probabilities Po and P1 before next time step both these indexes of are calculated by formula (1). First one scr directly for current situation and scr+ for virtual one when L[i,j] has one integer unit more.

Other method for reducing number of operations in algorithm is based on the optimization of updating procedure for matrices of probabilities to move of elements of matrix L, Po and P1. A list of all elements are in need to be updated, these includes all elements these changed position plus their neighbors, is created by algorithm at any time step. If size of the list is larger than number of all elements in active field of simulation all elements will be updated. Otherwise only elements in the list needed recalculation of probabilities are updated. This method produces most saving of calculation time for situation of to two dimensional nucleation, exactly when it is mostly important to make many time steps to achieve stable results.

The algorithm described above works close to theoretically achievable rate retaining complete integrity of underlined scientific formulas. The main acceleration compare to naive calculation with float values using exponent functions are coming because two major features: first one that it operates only with integer numbers and second that all repeated operations namely calculations of probabilities, these include exponential functions, are done in advance at preparation to simulation run with results stored in operative memory for fast reference as values of integer arrays.

The software algorithm of numerical simulation experiment for described above model also can be also summarized as sequence of following steps:

In spite of relatively straightforwardness of scientific background of the model the software implementation demands a few states of art solutions to increase effectiveness to the level of practical applicability. The critical point is to pre-calculated and save in operative memory most of intermediate values that permits saving of processor time from repeatedly calculated time consuming procedures needed for exponent functions, multiplication and division to use instead only integer sum and subtraction operations with indexes and values of integer matrices .

In details the algorithm can be learned through step-by-step examination its specific software implementation done in C++ codes project that can be provided on the base of case by case negotiations.

Oct. 24, 2017; 11:52 EST