A New Model to Predict Optimum Conditions for Growth of 2D Materials on a Substrate

Deposition of atoms or molecules on a solid surface is a flexible way to prepare various novel two-dimensional materials if the growth conditions, such as suitable surface and optimum temperature, could be predicted theoretically. However, prediction challenges modern theory of material design because the free energy criteria can hardly be applied to this issue due to the long-standing problem in statistical physics of the calculations of the free energy. Herein, we present an approach to the problem by the demonstrations of graphene and γ-graphyne on the surface of copper crystal, as well as silicene on a silver substrate. Compared with previous state-of-the-art algorithms for calculations of the free energy, our approach is capable of achieving computational precisions at least 10-times higher, which was confirmed by molecular dynamics simulations, and working at least four orders of magnitude faster, which enables us to obtain free energy based on ab initio calculations of the interaction potential instead of the empirical one. The approach was applied to predict the optimum conditions for silicene growth on different surfaces of solid silver based on density functional theory, and the results are in good agreement with previous experimental observations.


Introduction
Since graphene was obtained in 2004 [1], 2-dimensional (2D) materials have developed a wide interest all over the world. Due to the ultrathin thickness, the materials exhibit unique electrical, mechanical, and thermal properties, which inspires explorations of other more 2D materials. So far, dozens of 2D materials have been prepared experimentally, including the graphene family (e.g., graphdiyne [2], silicene [3], germanene [4], borophene [5], phosphorene [6], bismuthene [7]), III-nitrides [8,9], III-bismides [10], transition metal dichalcogenides [11], metal carbides [12], and the like, yet it remains a problem to prepare high-quality 2D material with larger sizes. Although a few kinds of 2D materials such as graphene and MoS 2 [13] can be obtained by exfoliating corresponding layered bulk materials, others, such as silicene [3,[14][15][16][17][18][19][20][21][22], cannot be produced in this way because there exists no corresponding layered bulk material in nature. Vapor deposition of atoms on a substrate should be much more flexible and has been applied for preparing various 2D materials of large scale [3,23]. However, much time and effort have to be paid to explore the growth conditions because the surface structure and the temperature of the substrate both have significant effects on the growth of 2D materials. As an example for preparing silicene by deposition of Si atoms on a silver substrate, the growth on the (001) surface only produces a "complex" superstructure without clear symmetry [24], and on the (110) surface, silicon nanoribbons (NRs) form along the [110] direction with a honeycomb structure [25,26], while on the (111) surface, a continuous graphene-like 2D honeycomb arrangement

Theoretical Method
The model for a 2D materials on a substrate is shown in Figure 1, where the substrate of M atoms is treated as a thermal bath at temperature T. For calculations of the PF for the 2D materials of N atoms, the total potential is expressed as: where U 2D is the potential energy of the 2D material with the coordinates of the atoms denoted by x 3N (x 1 , x 2 . . . x 3N ) and V is the interaction potential between the 2D material and the substrate with its atoms denoted by X 3M (X 1 , X 2 . . . X 3M ).
where 2 is the potential energy of the 2D material with the coordinates of the atoms denoted by 3 ( 1 , 2 … 3 ) and V is the interaction potential between the 2D material and the substrate with its atoms denoted by 3 ( 1 , 2 … 3 ).
The PF of the canonical ensemble for 2D materials can be expressed as: where β = 1/ with kB the Boltzmann factor and Q is the configurational integral: In order to solve the 3N-fold Q integral, the sense of the integral is reinterpreted as follows [35].
Traditionally, a 1D integral 1 = ∫ ( ) is interpreted as the sum of an infinite number of rectangles with area = ( )∆ , i.e., 1 = lim ∆ →0 ∑ . From another angle, the length of the 1D element ∆ at is modulated by ( ) to be a new length element ∆ ′ = ( )∆ and 1 = ∑ ∆ ′ . In other words, the 1D integral is a summation of length elements instead of area elements and equals an effective length of |b− |. Similarly, a 2D integral 2 = ∫ ∫ 0 0 ( , ) equals an effective area of • because the area element ds = dxdy is enlarged (or shrunk) by f(x,y), giving rise to an effective area element d ′ = ( , ) . Followed by this notion, an N-  The PF of the canonical ensemble for 2D materials can be expressed as: where β = 1/k B T with k B the Boltzmann factor and Q is the configurational integral: In order to solve the 3N-fold Q integral, the sense of the integral is reinterpreted as follows [35].
Traditionally, a 1D integral I 1D = b a f (x)dx is interpreted as the sum of an infinite number of rectangles with area A i = f (x i )∆x, i.e., I 1D = lim ∆x→0 i A i . From another angle, the length of the 1D element ∆x at x i is modulated by f (x i ) to be a new length element ∆x i = f (x i )∆x and I 1D = i ∆x i . In other words, the 1D integral is a summation of length elements instead of area elements and equals an effective length of |b−a|.
Similarly, a 2D integral I 2D = a 0 b 0 dxdy f (x, y) equals an effective area of a·b because the area element ds = dxdy is enlarged (or shrunk) by f (x,y), giving rise to an effective area element ds = f (x, y)dxdy. Followed by this notion, an N-fold integral I ND = When the integrand f (x 1 , . . x N ) being positive definite within the entire integral domain and having the minimum at the origin (U(0) = 0), the effective length of a i is defined as [35]: and the effective volume approximates to a product N i=1 a i , i.e., For the proof of Equation (5), we consider a 2D integral: where U(x, y) is positive definite within the integral domain, 0 ≤ x ≤ a, 0 ≤ y ≤ b, on which we define a map of x and y: Therefore, X (0) = 0, Y (0) = 0. The effective length of a and b is defined by: Inserting Equation (7) into Equation (6) yields: with: which can be expanded in Taylor series as: If U(x, y) has a minimum at the origin (x = 0, y = 0), i.e., U(0, 0) = 0, then F(0, 0) = 0, ∂F ∂X | 0 = ∂F ∂Y | 0 = 0, and thus, the value of F(X , Y ) is close to zero in the neighborhood of (0, 0). Since U(x, y) is positive definite, a and b are small, and the integral domain of Equation (9) is a small area around the origin, we obtain that: The proof can be easily extended to an N-dimensional integral, , as long as the function U(q 1 , q 2 . . . q N ) is positive definite and has a minimum at the origin (U(0) = 0). For the 3N-fold integral of Equation (3), although the integrand is of the same form as required by Equation (5), it may not be positive definite or have no minimum at the origin. Letting q 3N = q 1 , q 2 . . . q 3N be the coordinates of particles in the state of the lowest potential energy U 0 , we introduce a function: with x i = x i − q i . By inserting Equation (13) into Equation (3), we obtain: Clearly, U (x 3N , X 3M ) is positive definite within all the integral domains and has the minimum at the origin (U (0, X 3M ) = 0). According to Equation (5), the integral in Equation (14) equals an effective 3N-dimensional volume, where the effective length L i on the i th degree of freedom is defined as: For a 2D material sheet on a substrate ( Figure 1), however, the effective length L z might be different from L y or L x , and the edge atoms (N 1 ) should have different effective lengths from the ones of the atoms (N 2 ) in the center region. In such a case, Equation (15) turns into: To obtain the effective lengths, the first step is to find the most stable structure of the 2D materials with the lowest potential U 0 , which can be accomplished by various well-developed global optimization algorithms [39][40][41][42] or the dynamic damping method [43,44]. Starting from the most stable structure, one atom in the center region (or in the edge region) is moved step by step in one of the degrees of freedom, such as the X-axis, while the Y-and Z-coordinates and all other atoms are kept fixed to determine Clearly, it is an easy task for traditional ab initio algorithms and recently-developed DFT [45,46].
The computational cost of DIA and ANS [33] is determined by the number of times of the calculations of the potential energy. For a system consisting of hundreds of atoms [33,37,47,48], ANS partitions configuration space into at least 10 3 subdivisions, and in each subdivision, more than 3 × 10 4 configurations should be randomly produced to be calculated for the total potential energy, and the same program of ANS must be repeatedly run dozens of times to produce average results because of the fluctuations of the Monte Carlo algorithm in ANS. Therefore, the times for ANS calculating the potential are more than 3 × 10 8 , while the times needed by DIA can be fewer than 1 × 10 4 because (9) can be well determined with a step of 0.001 Å for changes of x i in a range of about 1 Å, indicating that DIA works at least four orders of magnitude faster than ANS.

Graphene and γ-Graphyne on a Cu Substrate
The approach developed in the last section was first demonstrated by calculating the PF for a piece of graphene (Figure 2a) or γ-graphyne ( Figure 2c) sheet of 510 C atoms on the (111) surface of a Cu substrate of 2640 atoms located in perfect fcc lattices. Brenner potential was employed to characterize C atoms' interaction with the empirical parameters taken from [34], and V(x 3N , X 3M ) was taken as the summation of the Lennard-Jones (L-J) potential f (r) = 4ε( σ 12 r 12 − σ 6 r 6 ) for the pairwise interaction between a C and a Cu atom with ε = 0.0168 eV, σ = 2.2 Å [49]. As shown in Figure 3, the internal energy ( ) derived from the PF was in excellent agreement with that ( ) obtained from the MD simulations with a relative standard deviation of 0.00002%, which is too small to be shown in Figure 3. In the temperature range from 100 K to 1300 K, the relative error ( ) as below 0.03%, and only 0.0005% and 0.002% for graphene at 500 K and γ-graphyne at 1100 K, respectively. We cannot make the comparisons of the precisions between DIA and that of ANS [33] because the implementation of ANS costs too many computer hours for this system and the followed ones in this work. Nevertheless, the internal energy errors for ANS applied in crystal argon systems of 500 atoms characterized by the L-J potential function, which is much simpler than the many-body Brenner-function used here, were above 10% [37], which was much larger than the above ones. It is notable that the dependence of internal energy (E) on temperature is nearly linear, indicating that = 0 + , where B is a constant. According to the statistical physics, the constant B equals three for 3D crystal atoms with harmonic coordinates. However, for the graphene and γ-graphyne sheet, the constant B equals 2.97 and 2.90, respectively, showing that the C atoms are not in the motion of harmonics.  (111) substrate. The potential energy U (b,d) felt by a C atom moving along the X-, Y-, or Z-axis depends on the specific surrounding of the C atom.
The system was cooled below 0.01 K by a damping method [43,44] to determine the lowest energy U 0 and the most stable structure. According to the configuration (Figure 2a,c), the C atoms were divided into a center group and an edge group, and for each of the atoms with different surroundings, one of its coordinates (such as x ) was changed step by step with an interval of 0.001 Å, while its y and z coordinates and all other atoms were kept fixed to record U (0 . . . x i . . . 0, X 3M ), as shown in Figure 2b and d for graphene and γ-graphyne, respectively. U for the center atom is indeed different from that of the edge atoms, and for the same atom, U along one coordinate axis is also different from that along the other two, so the configuration integral was calculated as: By applying Equation (2) and E = − ∂ ∂β ln Z, the internal energy (E PF ) was obtained through: with temperature difference ∆T = 0.1 K. In order to test the accuracy, a common procedure for MD simulations of a canonical ensemble [43,44] was employed to produce the internal energy of the 2D materials in contact with a Nanomaterials 2019, 9, 978 7 of 14 thermal bath at a given temperature T. Specifically, the atoms of the substrate ( Figure 1) were fixed while the C atoms moved according to classical equations of motion, which were solved by the Verlet algorithm with a time step of 0.2 fs. Within the first 400 fs, all the carbon atoms were assigned velocities every 40 fs according to the Maxwell velocity distribution at temperature T, and then, the internal energy (E MD ) and the temperature were recorded every 30 fs to perform the average over 100 records.
As shown in Figure 3, the internal energy (E PF ) derived from the PF was in excellent agreement with that (E MD ) obtained from the MD simulations with a relative standard deviation of 0.00002%, which is too small to be shown in Figure 3. In the temperature range from 100 K to 1300 K, the relative error ( × 100%) as below 0.03%, and only 0.0005% and 0.002% for graphene at 500 K and γ-graphyne at 1100 K, respectively. We cannot make the comparisons of the precisions between DIA and that of ANS [33] because the implementation of ANS costs too many computer hours for this system and the followed ones in this work. Nevertheless, the internal energy errors for ANS applied in crystal argon systems of 500 atoms characterized by the L-J potential function, which is much simpler than the many-body Brenner-function used here, were above 10% [37], which was much larger than the above ones.  Applying = − T ln , the free energy of graphene and γ-graphyne on Cu (111) was calculated (Figure 4a), showing that the free energy of graphene is always smaller than that of γ-graphyne in the temperature range from 100 K to 3000 K. The difference at 100 K, 321 eV, decreased gradually down to 267 eV with the temperature increasing up to 3000 K, indicating that graphene should be more easily grown than γ-graphyne on a Cu (111) surface via depositing C atoms. We also calculated the free energy on a Ni (111) surface (Figure 4b) using the same method and found that graphene still owned smaller free energy than γ-graphyne, although the difference at 100 K, 318 eV, became gradually small with the temperature. These results are consistent with previous experimental observations [50,51]. It is notable that the dependence of internal energy (E) on temperature is nearly linear, indicating that E = U 0 + BNk B T, where B is a constant. According to the statistical physics, the constant B equals three for 3D crystal atoms with harmonic coordinates. However, for the graphene and γ-graphyne sheet, the constant B equals 2.97 and 2.90, respectively, showing that the C atoms are not in the motion of harmonics.
Applying F = −k B T ln Z, the free energy of graphene and γ-graphyne on Cu (111) was calculated (Figure 4a), showing that the free energy of graphene is always smaller than that of γ-graphyne in the temperature range from 100 K to 3000 K. The difference at 100 K, 321 eV, decreased gradually down to 267 eV with the temperature increasing up to 3000 K, indicating that graphene should be more easily grown than γ-graphyne on a Cu (111) surface via depositing C atoms. We also calculated the free energy on a Ni (111) surface (Figure 4b) using the same method and found that graphene still owned smaller free energy than γ-graphyne, although the difference at 100 K, 318 eV, became gradually small with the temperature. These results are consistent with previous experimental observations [50,51].
γ-graphyne in the temperature range from 100 K to 3000 K. The difference at 100 K, 321 eV, decreased gradually down to 267 eV with the temperature increasing up to 3000 K, indicating that graphene should be more easily grown than γ-graphyne on a Cu (111) surface via depositing C atoms. We also calculated the free energy on a Ni (111) surface (Figure 4b) using the same method and found that graphene still owned smaller free energy than γ-graphyne, although the difference at 100 K, 318 eV, became gradually small with the temperature. These results are consistent with previous experimental observations [50,51].

Silicene on Silver Substrate
For a piece of silicene sheet of 336 Si atoms on the (111) surface of an Ag substrate with 2640 Ag atoms located in perfect fcc lattices (Figure 5a), the Tersoff potential was employed to describe the interactions between Si atoms with the parameters taken from [36], and the interaction between the

Silicene on Silver Substrate
For a piece of silicene sheet of 336 Si atoms on the (111) surface of an Ag substrate with 2640 Ag atoms located in perfect fcc lattices (Figure 5a), the Tersoff potential was employed to describe the interactions between Si atoms with the parameters taken from [36], and the interaction between the Si atoms and the Ag atoms as described by a Morse pairwise function [52]. The system was cooled down to 0.01 K by a damping method [43,44] to determine the lowest energy U 0 and the most stable structure. Although U (0 . . . x i . . . 0, X 3M ) felt by an atom i located at the edges of the silicene sheet should be different from the one felt by an atom in the center region, the edge atoms are much fewer than the center atoms, so we only calculated U felt by a center atom (Figure 5b) to determine the effective length L x , L y , and L z , and the configuration integral approximates to: according to Equation (2) and E = − ∂ ∂β ln Z, the internal energy (E PF ) can be obtained through Equation (19) with ∆T = 0.1 K.
A classical MD simulation was performed by using the Large-scale Atomic/Molecular Massively Parallel Simulator [53] program to produce the internal energy of the system with the same potentials [36,52] as those for calculations of the PF. Specifically, the time step as set as 0.1 fs, and the velocities of the silicon atoms were assigned according to a time integration on Nose-Hoover equations of motion to keep the system at a given temperature. Then, the internal energies (E MD ) and the temperature were recorded every 30 fs to perform the average over 100 records.
As shown in Figure 6, the internal energy (E PF ) derived from the PF was in excellent agreement with that (E MD ) obtained by the MD simulations with a relative standard deviation of 0.008%, which is too small to be shown in Figure 6. For temperatures lower than 500 K, the relative error ( ) was only 0.001% and gradually increased up to 0.021% for 1300 K. It may be expected that the relative error would get smaller if the effective length L i of the edge atoms was calculated, instead of replacing L i with the one of the center atom, to calculate Q by Equation (11).
A classical MD simulation was performed by using the Large-scale Atomic/Molecular Massively Parallel Simulator [53] program to produce the internal energy of the system with the same potentials [36,52] as those for calculations of the PF. Specifically, the time step as set as 0.1 fs, and the velocities of the silicon atoms were assigned according to a time integration on Nose-Hoover equations of motion to keep the system at a given temperature. Then, the internal energies ( ) and the temperature were recorded every 30 fs to perform the average over 100 records. As shown in Figure 6, the internal energy ( ) derived from the PF was in excellent agreement with that ( ) obtained by the MD simulations with a relative standard deviation of 0.008%, which is too small to be shown in Figure 6. For temperatures lower than 500 K, the relative error (

Conditions for Silicene Growth on a Ag Substrate
In order to search for the optimum conditions for given silicene grain growth on a Ag substrate, we calculated the FE ( = − ) by DIA with the potential ′ determined by DFT, which were performed by the Vienna ab initio Simulation Package based on local density approximation. The kinetic-energy cutoff for the plane-wave basis set was 400 eV, and the Brillouin zone was sampled with (2 × 2 × 1) k-points.
Considering that silicene can grow on a Ag (110) surface by deposition of Si atoms, we calculated the FE of a silicene grain consisting of 35 Si atoms in hexagon arrangement with one Si atom adhered to the zigzag (Figure 7a) or armchair edge (Figure 7b) on the surface of four atomic layers with each containing 40 Ag arranged in an 8 × 5 fcc supercell. According to thermodynamics, if the free energy for the zigzag adherence (FEZ) equals the one for the armchair adherence (FEA), then both the zigzag and armchair edge of the grain are able to grow larger. Otherwise, the future Si atoms deposited on the Ag surface would more favorably arrive at the zigzag (or armchair) edge if FEZ is lower (or larger)

Conditions for Silicene Growth on a Ag Substrate
In order to search for the optimum conditions for given silicene grain growth on a Ag substrate, we calculated the FE (F = −k B T ln Z) by DIA with the potential U determined by DFT, which were performed by the Vienna ab initio Simulation Package based on local density approximation. The kinetic-energy cutoff for the plane-wave basis set was 400 eV, and the Brillouin zone was sampled with (2 × 2 × 1) k-points.
Considering that silicene can grow on a Ag (110) surface by deposition of Si atoms, we calculated the FE of a silicene grain consisting of 35 Si atoms in hexagon arrangement with one Si atom adhered to the zigzag (Figure 7a) or armchair edge (Figure 7b) on the surface of four atomic layers with each containing 40 Ag arranged in an 8 × 5 fcc supercell. According to thermodynamics, if the free energy for the zigzag adherence (FEZ) equals the one for the armchair adherence (FEA), then both the zigzag and armchair edge of the grain are able to grow larger. Otherwise, the future Si atoms deposited on the Ag surface would more favorably arrive at the zigzag (or armchair) edge if FEZ is lower (or larger) than FEA, resulting in development of only one of the edges to form a band (or nanoribbon). For calculations of FEZ and FEA, the system was optimized firstly with the bottom layer of the Ag substrate fixed, and then, the adhered Si atom (and one of the Si atoms denoted by a square in Figure 7a and b in the center of the silicene) was moved along the X, Y, and Z direction step by step with an interval of 0.1 Å to produce ′ for calculating the effective length , , and (or the , , and ) via Equation (9), and finally, the configurational integral was obtained by: where 0 = 35, the number of Si atoms except for the adhered Si atom. As shown in Figure 7c, FEZ and FEA decreased significantly with the temperature increasing up to 2000 K, while the difference DFE (= FEZ-FEA) decreased gradually from 2.854 eV for 0 K down to 2.526 eV for 2000 K. According to the thermodynamics, if a Si atom deposited on this surface has enough time to wander between the two edges, then it will finally locate at the armchair edge because of the lower FE, and five similar Si atoms will form an armchair edge instead of a zigzag one. As a result, the armchair edge progresses row by row, and the initial silicene grain eventually grows into a nanoribbon along the [110] direction with a width of about 1.6 nm, which is the exact observation in the previous experiment [25,26].
For the growth on the Ag (111) surface, a silicene grain of 25 Si atoms arranged in hexagons with a Si atom adhered at the zigzag edge (Figure 8a) or the armchair edge (Figure 8b) was placed on an Ag substrate of four (111) layers with each consisting of 48 Ag atoms in a (6 × 8) supercell, and the geometry was optimized with the bottom layer of the Ag substrate fixed. The adhered Si atom (and the Si atom denoted by a square in Figure 8a,b) was moved step by step to produce the potential ′ , and the configurational integral was obtained by Equation (21) with 0 = 25. As shown in Figure 8c, the FEZ and FEA decreased with the temperature up to 2000 K, while the DFE increased slightly from −0.022 eV up to 0.056 eV. In the temperature range from 200 K to 580 K, the DFE as less than 0.01 eV, and the DFE equaled zero for 400 K. According to the thermal dynamics, the probability for the future deposited atoms adhering to either the zigzag or the armchair edge was equal when the substrate was kept at 400 K, which should be the most optimum temperature for the silicene grain growing to form a continuous graphene-like structure. In previous experiments [3,27], the optimum temperature was 493-580 K, which was about one hundred Kelvin higher than our prediction. The difference may have resulted from the limited computational precision of DFT. For calculating 0 and ′ (0 … ′ … 0, 3 ) in Equation (13) and (16), it is well known that the calculation results of DFT For calculations of FEZ and FEA, the system was optimized firstly with the bottom layer of the Ag substrate fixed, and then, the adhered Si atom (and one of the Si atoms denoted by a square in Figure 7a,b in the center of the silicene) was moved along the X, Y, and Z direction step by step with an interval of 0.1 Å to produce U for calculating the effective length l x , l y , and l z (or the L x , L y , and L z ) via Equation (9), and finally, the configurational integral was obtained by: where N 0 = 35, the number of Si atoms except for the adhered Si atom. As shown in Figure 7c, FEZ and FEA decreased significantly with the temperature increasing up to 2000 K, while the difference DFE (= FEZ-FEA) decreased gradually from 2.854 eV for 0 K down to 2.526 eV for 2000 K. According to the thermodynamics, if a Si atom deposited on this surface has enough time to wander between the two edges, then it will finally locate at the armchair edge because of the lower FE, and five similar Si atoms will form an armchair edge instead of a zigzag one. As a result, the armchair edge progresses row by row, and the initial silicene grain eventually grows into a nanoribbon along the [110] direction with a width of about 1.6 nm, which is the exact observation in the previous experiment [25,26].
For the growth on the Ag (111) surface, a silicene grain of 25 Si atoms arranged in hexagons with a Si atom adhered at the zigzag edge (Figure 8a) or the armchair edge (Figure 8b) was placed on an Ag substrate of four (111) layers with each consisting of 48 Ag atoms in a (6 × 8) supercell, and the geometry was optimized with the bottom layer of the Ag substrate fixed. The adhered Si atom (and the Si atom denoted by a square in Figure 8a,b) was moved step by step to produce the potential U , and the configurational integral was obtained by Equation (21) with N 0 = 25. As shown in Figure 8c, the FEZ and FEA decreased with the temperature up to 2000 K, while the DFE increased slightly from −0.022 eV up to 0.056 eV. In the temperature range from 200 K to 580 K, the DFE as less than 0.01 eV, and the DFE equaled zero for 400 K. According to the thermal dynamics, the probability for the future deposited atoms adhering to either the zigzag or the armchair edge was equal when the substrate was kept at 400 K, which should be the most optimum temperature for the silicene grain growing to form a continuous graphene-like structure. In previous experiments [3,27], the optimum temperature was 493-580 K, which was about one hundred Kelvin higher than our prediction. The difference may have resulted from the limited computational precision of DFT. For calculating U 0 and U (0 . . . x i . . . 0, X 3M ) in Equations (13) and (16), it is well known that the calculation results of DFT depend significantly on the specifically employed basis sets and exchange-correlation functionals of the electron density, for which the recent work [45,46] might provide better choices.

Summary
In summary, DIA was developed to calculate the free energy (or PF) of 2D materials on a substrate, and the high calculation precision was validated by MD simulations. It should be pointed out that such a test is much more stringent than the comparisons between the results derived from the PF and experiments because the same interaction potential is used in both calculations of the PF and the MD simulations, while the potential may not correctly describe the realistic interactions between atoms concerned in the experiment. As for the efficiency, DIA works at least four orders of magnitude faster than the most efficient method, ANS, developed previously, and enables calculations of the free energy based on ab initio calculations to predict the optimal conditions for novel 2D materials' growth on a substrate, which would greatly shorten the way to experimental realization.

Summary
In summary, DIA was developed to calculate the free energy (or PF) of 2D materials on a substrate, and the high calculation precision was validated by MD simulations. It should be pointed out that such a test is much more stringent than the comparisons between the results derived from the PF and experiments because the same interaction potential is used in both calculations of the PF and the MD simulations, while the potential may not correctly describe the realistic interactions between atoms concerned in the experiment. As for the efficiency, DIA works at least four orders of magnitude faster than the most efficient method, ANS, developed previously, and enables calculations of the free energy based on ab initio calculations to predict the optimal conditions for novel 2D materials' growth on a substrate, which would greatly shorten the way to experimental realization.