Physical background of the numerical model for crystal growth is described in details. There are definition of the simulation field and formulas for probabilities of elemental moves of molecules.
A model of growing crystal has to incorporate adequate mathematical presentation of the crystal and molecule moves applicable for software implementation. It has to reflect the molecule structure of the crystal surface including symmetrical positioning of molecules together with rules derived from atomic physics for calculation probabilities for each molecule to change their positions in small time step. Following is a description of the model used in software application LeoMonteCrystal.
Simulation field.
The physical reality of growing orthorhombic crystal is represented in the model by an integer matrix L[N_{b}, N_{c}] whose values are positions of the molecules on crystal surface measured in layers of crystallographic lattice in direction of the growth. Indexes of the matrix L give location of the projection of a molecule on the simulation field, where its dimensions, N_{b} and N_{c}, give a size of the simulation field. Fig.1 shows how for each element of matrix L these three numbers, one of them is for value of matrix and two others are for its indexes, unambiguously define a position of the molecule on the surface of the crystal. A value of the matrix element provides distance from base layer that also can be chosen as a starting position for the numerical experiment of growing surface. Two indexes represent projection of molecule on the base layer in the Cartesian coordinate system oriented consistently with orthorhombic lattice of the crystal. The positions of molecules in physical space can be obtained from integer values of matrix and its integer indexes, when all three of them are measured in layers of crystal, by scaling them by corresponding constants of crystallographic lattice.




Fig. 1. Visual presentation how two dimensional integer matrix represents a crystal surface in three dimensional space. Borders of active simulation field of matrix L are highlighted by blue color. 

To avoid complications dealing with special conditions on edges of matrix L, elements at opposite sides of active field are glued to each other creating infinite kaleidoscope like space of simulation. It makes element L[i,N_{c1}] neighbor to element L[i,0] and L[N_{b1},i] is touching L[0,j] (here and later we are using matrix L when zero based index notation is applied as it customs for C++ language that is used for software implementation).
In software implementation of the algorithm matrix L has two plus rows and columns: L[N_{b}+2, N_{c}+2] so next state of matrix L at first calculated on the active part of the matrix that doesn't include border rows and columns. Elements of these marginal rows and columns are taking into consideration nevertheless for calculating neighborhood of molecules inside of active field. After defining new values of elements in the active field these border elements are assigned to corresponding elements at active simulation area by following rules:
for margin columns:  L[0,i]=L[N_{b},i] and L[N_{b}+1,i]=L[1,i]; 
for margin rows:  L[j,0]=L[j,N_{c}] and L[j,N_{c+1}]=L[j,1] 
In the Fig.1 borders of active simulation field of matrix L is highlighted by blue color.
Characteristics of the crystal surface these are in correspondence to the state of the matrix L can be calculated by the following formulas:
Position of the surface is defined as average value of all elements matrix L belonged to active field. It can be calculated by formula:
(1) 
where summing is performed for all n = N_{b} *N_{c} elementsof active field of the matrix L.
Average growth rate is defined by formula:
(2) 
where t_{1} and t_{2} are time checkpoints (arbitrary chosen) to calculate average rate of growth in between two time points during simulation experiment.
Roughness of the growing surface is calculated as standard deviation position of molecules from the average position P:
(3) 
If all elements in matrix L have the same value it means that the surface is 100% flat (absolute smooth). For absolutely flat surface the value of roughness equals zero (S =0.0). There is no theoretical limitation on the upper side of of the roughness parameter except applying general physical sense that suggests that more than ten layers roughness could be unrealistic.
Concentration of hole like defects below of the growing surface is calculated as ratio of number evens of creating hole N_{h+} minus event to open existing hole N_{h} to all position above base layer P*n.:
C_{h} = (N_{h+ }  N_{h})/(P*n)  (4) 
The elementary event of creation hole like defect is formally associated with act of accepting a molecule at special bottle neck (well) like position in the matrix L. Such bottleneck position in matrix L_{i} is formally defined as having value by at least two less than all its closest neighbors. There is presumption that when molecule is incorporated in such position it does not move up to the bottom of "well" and leaves empty space below that (empty space) is interpreted as hole like defect. There is also an option of molecule above hole to be emitted from crystal open the hole is opened and a number of defects is lessened by one. During simulation every such elementary act is arithmetically summed to got actual number of hole defects below growing surface of the crystal.
An algorithm for numerical simulation is based on Monte Carlo implementation of modified cellular automaton model. The state of each element of matrix L for the next time step is changed randomly proportionally to calculated probabilities for corresponding molecule to make a jump into crystal or from it. These probabilities are calculated separately for each position at matrix L and is depending on relative positions of near neighbor molecules. For each molecule its nearest neighbors surrounding defines effective change of the surface of crystal in case of emitting molecule from crystalline phase into amorphous (noncrystalline) feeding matter surrounded crystal. A simulation algorithm can be described as infinite loop of calculations of next state of matrix L based on its current situation by set of rules these define probabilities for any given element of the matrix after certain period of time to change its position. There are three options for molecule at any give position after given time step, with the value set by researcher, to move ahead, back or stay at the same position. By application of random numbers generator one of these options is chosen proportionally corresponding probabilities for each element of the matrix to change position. That procedure will define state of matrix L at the next time step. Then this procedure repeats itself to calculate state of matrix L for next time step in infinite loop until reliable statistics of the process should be obtained.
Each element of matrix L represents corresponding position the surface of the crystal. At any given time step the value of element has an a choice to:
First three options have probabilities p_{+}, p_{} and p_{o} these are different for each element of matrix L. An actualization of the option number four (hole like defect) can be checked after accepting next molecule or emitting one. At any given time step in algorithm there are also two float (or in fact normalized to maximum random number integer) matrixes with the same size as matrix L for following probabilities:
Calculation of probabilities for move of each element.
In accordance to the theory of thermally activated reactions the frequency of incorporation of a molecule into crystal:
γ_{+} = γ_{o} • _{e}^{(E a in/kT) }  (5) 
and an opposite movement of molecule from the crystal into noncrystalline phase:
γ_{} = γ_{o} • _{e}^{(E a out/kT) }  (6) 
are defined by modified Arrhenius equation where γ_{o}  a frequency of thermal vibrations of a molecule, T  absolute (Kelvin scale) temperature, k  Boltzmann's constant (1.380 65×10^{23} J K^{1}), E_{a in }and E_{a out }– energies for activation barriers for process of incorporation of the molecule into crystal and for opposite process of emitting it from crystal body into feeding phase correspondingly. The theory of thermally activated reactions postulates that molecule on surface to be incorporated into crystal must on the way to be deformed in some manner rearranging atoms into unstable intermediate atomic complex. Formation of such intermediate atomic complex (sometimes named activating complex) is accompanied by increase of free energy level during the jump (transformation) associated with the molecule. It's worth to note that jumps from crystalline phase and back not necessary associated with comparable to crystallographic lattice changing location of molecule in space on the scale of size of crystal layer. It can be accomplished by rearranging positions and bonds for neighboring atoms in location of molecule instead. The minimum necessary deformation is defined here is as such that demands a minimum of increasing energy level that permits to fulfill a phase transformation for single molecule.
A difference between energy levels of molecules in crystal and feeding phase is:
E_{a out } E_{a in } = ΔG•v_{m}  σ•Δs_{m}  (7) 
where ΔG is a change of specific Gibbs potential of free energy defined by equation:
ΔG = ΔH* (ΔT / T_{l})  (8) 
where ΔH – specific enthalpy of crystallization, ΔT
= T_{l} – T, where T_{l }– an equilibrium temperature of
the phase transition, T  temperature,
σ  surface energy of the crystal surface in contact with
surrounding phase, Δs_{m}  change of surface of
the crystal after emitting of one molecule from any given position
on the crystal (note that for absolute rough surface in case
of continuous mechanism this value should be zero at average
and for smooth surface at layer by layer growth value of Δs_{m
}should be significant large), v_{m} – volume of molecule.
The implementation of the discussed model demands establishing a connection between frequencies of jumping of molecules in and out of crystal: γ_{+} , γ_{} with probabilities for molecule on surface in given position of matrix L to move ahead, back or stay on the same position: p_{+}, p_{} and p_{o}.
Time step.
Average number of molecules that will be incorporated into surface of the crystal absent of the phenomena of emitting back, n_{+} , can be arbitrary defined by researcher as a starting point in defining value of time step. The value n_{+} should be small enough for adequate modeling of the process that should be achieved if random positions of these molecules n_{+ }will fall far enough from each other not to be touching neighbors for every step. Other way around, there is the need for time step to be large enough to fit into always limiting computing resources. For specificity sake one can consider value about n_{+} = 10 on the simulation field with size n = 100*100 as quite proper compromise for both these demands for modern computers. The physical value of time step can be defined from the value of frequency of the processes of incorporation of molecule into crystal by formula:
τ = (n_{+}/n)/μ_{+}  (9) 
Probabilities for surface changing position on the molecular scale.
If value of time step approaches zero, τ →0, probabilities of given position on surface of the crystal to move on one step back or forth, or stay where it is can be calculated by equations:
p_{+} = γ_{+} • τ = (n_{+}/n)  (10) 
p_{} = γ_{} • τ = p_{+}•exp((ΔG•v_{m} + σ•s_{m}•Δs_{m})/kT) = p_{+}•exp(ΔG`  σ`•Δs_{m})  (11) 
p_{o} = 1  (p_{+} + p_{})  (12) 
where general supercooling is defined as:
ΔG` = ΔG•v_{m}/kT  (13) 
and specific surface energy:
σ` = σ•s_{m}/kT  (14) 
The relative change of the crystal surface for event of emitting one molecule gives the formula:
Δs_{m} = (1  2 *s_{cr})  (15) 
where s_{cr} is part of molecule on surface that is in direct contact with near neighbor molecules these are belonged to the crystalline phase.
Absent of limitations on computer's productivity and numerical precisions in calculations including random number generator formulas (10  12) can be used straightforwardly if small enough time step if chosen. Understanding that for Δs_{m} → s_{m} the value p_{}, calculated by formula 10, has largest possible value and arbitrary choosing p_{}(s_{m})_{ } ≡ 0.5 or lesser, we will be sure that any probabilities calculated by formulas (10  12) with
n_{+} ≡ 0.5*n*exp(ΔG`  σ`)  (16) 
will be lesser than this arbitrary number (0.5). All calculated probabilities for movements of the molecules will be less than 1.0 and so will satisfy all the rules. There is large set of situations with proper combination of values of general supercooling and specific surface energy when such direct approach is perfectly workable. But not for every situation.
There are also infinite and most interesting set of big values of specific surface energy when this naive, direct approach, described above, cannot be realistically implemented even for modern computers. The problem manifests itself clearly in the fact that there always could be situations for large values of σ` when n_{+} should be too small for computers to handle. If value n_{} for described direct approach had happened to be very small it produces at least too complications: first one is that time of numerical experiment will too consuming even most powerful computers to achieve practically applicable results. Second problem, that could even more problematic, occurs when precision of the generator of random numbers that can be defined as inverse value of its maximum produced number, 1/N_{r_max}, ought to be less or close to probability of one molecule jump into crystal per one time step. That condition can be presented by formula:
n_{+}/n ≥ 1/N_{r_max}  (17) 
We use two methods to expand area of workable ΔG`, σ` values beyond limitations presented by formulas 16 and 17.
A part of molecule belonged to crystalline phase for given element of matrix L.
From previously presented formulas it is clear conclusion that parameter s_{cr}, part of the molecule surface in direct contact with its neighbor molecules these are belonged to crystalline phase, is a primary parameter for calculating probabilities of changing position of any given element of matrix L at next time step.





Fig. 2. Directions to near neighbors of molecule incorporated into a cell of orthorhombic lattice. 

As it is shown in the Fig.2 a molecule inside of the cell of orthorhombic crystallography lattice has 26 unequally close neighbors along directions from its center toward its 6 faces, 8 vertices and 12 edges (6+ 8 + 12 =26). If surface of molecule are divided to areas associated with each neighbor molecules the value of molecule surface can be presented by equation:
S_{m} = 2*(S_{a} + S_{b} +S_{c}) +4*(S_{ab} + S_{ac} + S_{cb}) + 8 *S_{abc}  (24) 
where S_{a}, S_{b}, S_{c} – effective parts of molecule surface toward corresponding faces; S_{ab}, S_{ac}, S_{cb} edges; and S_{abc} vertices.
There is other way to define surface of molecule as scaled surface of incorporated into cuboid cell molecule:
Note that in specific implementation of our software LeoMonteCrystal user can define these parameters via two ratios toward of faces part of surface: S_{e}= s_{e}/s_{f} and S_{v}= s_{v}/s_{f} in form following ratios 1 : S_{e }: S_{v}. In reverse taking into account equation 27 parts all relative parts of surface in these direction can be derived from these two ratios as following:
s_{f } = 1/(1 +S_{e }+ S_{v})

(28) 
It is logical to presume that relative input of different faces in whole molecule surface is proportional their individual surface:
S_{a} : S_{b} : S_{c} = b*c :^{ }a*c : a*b  (29) 
In their turn we proclaim that relative surfaces associated with edges are proportional to their lengths:
S_{ab} : S_{ac} : S_{cb} = c :^{ } b : a  (30) 
All parts of surface associated with vertices have the same value.
Summary. Matrix L permits adequate representation of the surface of orthorhombic crystal. On the base the theory of thermally activated reaction probabilities of the molecule to be incorporated and be emitted from crystal can be calculated for small enough time step for each molecule for any given position of molecule associated with element of matrix L taking into account values of the closest neighbors of the element to obtain values of molecule's surface in direct contact with crystal phase. A model of quasiindependent behavior of molecules during time step permits to expand areas of thermodynamic parameters in numerical simulation applicable for series of experiments run on modern computers for reasonable time frame.
Oct. 19, 2017; 12:46 EST