A review of atomic layer deposition modelling and simulation methodologies: Density functional theory and molecular dynamics

The use of computational modelling and simulationmethodologies has grown in recent years as researchers try to understand the atomic layer deposition (ALD) process and create new microstructures and nanostructures. This review article explains and simplifies two simulation methodologies, molecular dynamics and the density functional theory (DFT), in solving atomic layer deposition problems computationally. We believe that these simulation methodologies are powerful tools that can be utilised in atomic layer deposition. DFT is used to solve problems in surface science and catalysis (predicting surface energy, adsorption energy, charge transfer, etc.), semiconductors (band structure, defect bands, band gap, etc.), superconductors (electron–phonon coupling, critical transition temperature), and molecular electronics (conductance, current–voltage characteristics). Molecular dynamics (MD) is used to predict the kinetic and thermodynamic properties of a material. Of interest in this article is a review where different material problems emanating from atomic layer deposition from these fields have been addressed by DFT and MD. Selected publications are discussed where DFT and MD have been successfully applied in atomic layer deposition (and related processes in some instances). The applications of DFT stretch from binding energy calculations of molecules and the solid band structure in chemistry and physics, respectively, computing the electron density up to determining the properties of a many-electron system. Also highlighted in this review study are the challenges that DFT and MD simulations must overcome.


Introduction
Atomic layer deposition (ALD) is a growing material deposition method, used particularly where precise thickness control of the material growth is required [1,2]. It can be used to produce ultrathin films [3,4]. ALD can be studied and analysed experimentally and by computational modelling and simulations. Experiments are expensive to run since the ALD equipment and other raw materials are expensive. They can also pose hazards because sometimes gases used or released in the chemical reactions are toxic and occasionally explosive. Modelling and simulations reduce drastically the costs involved in experiments. ALD [2,[5][6][7][8][9] has grown to become a vital manufacturing technique in many fields of nanotechnology, [10][11][12] semiconductor processing [13][14][15][16], microelectronics [17][18][19][20][21], electrocatalysis [22][23][24], and protective coatings [25][26][27][28][29][30][31]. ALD is a very powerful method of fabricating nanostructures and nanomaterials, regardless of their complex geometric shapes, with uniform films and coatings at the nano, micro, and atomic levels [7,13,25,[32][33][34]. The computer simulation of molecular systems has recently increased in academic and industrial research [35]. Computational modelling and simulation methodologies have several advantages over experiments. Some of these are low costs because they do not require costly equipment and raw materials. Because of their versatility, they may be used to investigate and analyse practically any material virtually, including those that produce toxic gases during their manufacture. They can be performed almost anywhere if adequate computational power is available. Atomic layer deposition can be modelled and simulated using various modelling and simulation methodologies such as computational fluid dynamics (CFD), density functional theory (DFT), Knudsen number, Monte Carlo method, and molecular dynamics (MD). ALD can be analysed at atomic, nano, and micro scales. It is essential to know which methodology one must adopt in modelling and simulating ALD. This review article looks at the two primary modelling and simulation methodologies that are employed in ALD: the DFT and molecular dynamics. Molecular dynamics uses Newton's laws of motion to simulate forces of attraction and the position of molecules, while the DFT uses quantum mechanics principles to simulate the interaction of molecules and atoms. The problems that one can expect to solve using DFT in atomic layer deposition are predicting the surface energy of a material, adsorption energy calculations, charge transfer predictions, band structure prediction in semiconductors, defect bands identification, band gap calculation in semiconductors, electron-phonon coupling prediction in superconductors, critical transition temperature prediction in superconductors, and conductance and current-voltage characteristics in molecular electronics. MD is used to predict the kinetic and thermodynamic properties of a material. Other review articles have discussed DFT/MD for ALD. These include the studies by Oviroh et al. [2], where they discuss various developments in ALD and why researchers opt for computational methodologies, Elliott et al. [36], where various classes of ALD processes based on mechanistic information from DFT are discussed, and Gupta et al. [37], who review various advances in materials design using ALD in the energy application fields. The review will also consider where the methods have been applied, their pros, and challenges to date.
2 Atomic layer deposition process ALD [3,5,25,[38][39][40][41][42][43][44] is a film deposition technique used to synthesise ultra-thin films. The process is based on sequential self-limiting chemical reactions. In an ALD process, a substrate on which thin films must be deposited is placed inside a reaction chamber at a specific pressure and temperature. The first precursor is then injected into the reaction chamber and reacts with the surface of the substrate forming a thin layer. By-products and excess first precursor molecules are eliminated by blowing in inert gas. The second precursor is then blown in. A chemical reaction occurs between the first precursor surface on the substrate and the second precursor. The excess precursor molecules and by-products are also carried away by the inert gas blown in after this reaction. The cycle of precursor 1 → purge → precursor 2 → purge → precursor 3 → purge, is repeated until the required film thickness has been achieved. Figure 1 shows the illustration of the ZnO atomic layer deposition using two precursors (C 2 H 5 ) 2 Zn (zinc source) and H 2 O (oxygen source).
There are several variations of atomic layer deposition, namely, metal ALD [21], thermal ALD [45,46], particle ALD [41], plasma-enhanced ALD [45,[47][48][49], photoassisted ALD, and area selective ALD [34,[50][51][52]. Metal ALD is used to deposit metals on a substrate by an elimination reaction. Most metals, except aluminium, can be deposited by metal ALD. In most cases, the precursors are a metal halogen (usually a fluoride) and a silicon-based precursor. Temperatures between 175 and 325°C are usually employed in metal atomic layer deposition. Thermal ALD is used where relatively high temperatures are required, for example, molecules that contain aluminium. Thermal ALD is used in Al 2 O 3 ALD using trimethylaluminium and water as precursors. The temperatures range between 150 and 350°C in thermal ALD. Particle ALD is closely related to the conventional atomic layer deposition process but is more concerned about coating the entire particle surface, including the surfaces of nanoparticles. It utilises the principles of ALD and applies them to the complex geometry coating of particles that are spherical in nature. The temperature range for this process varies with the nature of precursors and substrates but is usually lower than that of thermal and metal ALD. Plasma-enhanced ALD uses plasma to reduce the temperature required to vapourise molecules. Kim [53] reports that the following metals, metal oxides, and nitrides have been successfully deposited using PEALD: Ti, Ta, Al, Ru, Cu, Co, Ni, W, TiN, TaN, RuTiN, TiSi x N y , TaSi x N y , W 2 N, SiN x , Al 2 O 3 , Ta 2 O 5 , Y 2 O 3 , ZrO 2 , ZnO, HfO 2 , SrTiO 3 , SrTa 2 O 6 , and SrBi 2 Ta 2 O 9 . Photo-assisted ALD uses ultra-violet (UV) light to start and speed up surface reactions on the substrate. Photo-assisted ALD can be controlled by altering the wavelength, illumination time, and intensity of the ultra-violet light. Miikkulainen et al. [54] reported that few materials that include GaAs, BN, and metal oxides have been successfully deposited using photo-assisted ALD.
The ALD process has gained interest and appreciation among researchers in the past decade due to its precise and uniform atomic thickness control of the deposited film [2]. Among the different thin film fabrication processes, ALD is steadily growing to be the preferred thin-film deposition technique. Atomic layer deposition has several advantages over the other fabrication techniques, like precise atomic control of deposited film thickness, homogeneousness of thicknesses deposited even on complex three-dimensional surfaces, relatively lower system temperatures, and improved electronic properties of devices manufactured from films deposited by ALD [2]. Tables 1 and 2 compare the atomic layer deposition process with other thin-film fabrication techniques. Table 3 outlines the strengths and weaknesses of atomic layer deposition.

Atomic layer deposition in nanostructured materials
Nanostructured materials [2,20,32,33,59] are materials whose structural elements have dimensions that range from 1 to 100 nm [60]. These structural elements include molecules, crystallites, and clusters. Nanostructured materials can be categorised into four classes: carbon-based, inorganic-based, organic-based, and composite-based nanomaterials. Jeevanandam et al. [61] found that there are also other ways of classifying nanostructures such as their dimensions (0D, 1D, 2D, and 3D) and based on their origins (natural, which are produced in nature by anthropogenic activities or biological species, and synthetic/ engineered nanomaterials). In this review, we discuss briefly the ALD process on carbon-based and inorganicbased nanostructures.

ALD on carbon-based nanostructures
Classified as either metallic or semiconductors, they encompass many low-dimensional carbon allotropes like carbon black (CB), graphene, fullerene, carbon nanotubes (CNTs), and carbon fibre [62]. Carbon nanotubes are the strongest one-dimensional materials. Graphene has good electrical conductivity, high tensile strength, high transparency, and is the thinnest known two-dimensional material. Nanomaterials are applied in conducting materials, transducers, biochemical sensors, electronics, and conducting materials. ALD has been used to coat and decorate nanostructures and has proved to be the coating technique of choice for carbonbased nanostructures because of its conformal and self-limiting surface deposition. The first use of ALD to coat carbonbased nanostructures, in this case, carbon nanotubes, was reported by Javey et al. [63] in 2002. The group deposited high-k dielectric zirconium oxide thin films onto singlewalled carbon nanotubes. Lee et al. [64] reported the homogeneous coating of multi-walled CNTs (MWCNTs) with Al 2 O 3 from water and trimethylaluminium (TMA) precursors. Figure 2 shows the transmission electron microscope image of the Al 2 O 3 -coated MWCNTs.

ALD on inorganic-based nanostructures
Examples of inorganic-based nanostructured materials include transition metal oxides, carbon nanostructures, and metal nanoparticles. They comprise a non-metal or metal element and can also take the form of a hydroxide, oxide, chalcogenide, or phosphate compound [65,66].
They show some characteristics of electrical conductivity and are used in composite and hybrid materials for flexible electronic applications. Atomic layer deposition has been applied in manufacturing inorganic thin films using a 3D nanonetwork structure made of a polymer as a substrate. This substrate is usually pre-formed using advanced 3D nanofabrication techniques such as direct laser writing (DLW), lithography, block-copolymer (BCP), electrospinning, phase-mask interference lithography (PMIL), and multibeam interference lithography (MBIL). These 3D nanofabrication processes chiefly yield polymer-based nanonetwork structures. Before these polymers can be applied in various fields, they need to be improved. These polymers have poor mechanical resistance and poor thermal resistance. Typical methods used to fill the interconnected pores are electrochemical deposition, chemical vapour deposition, and sol-gel reaction. Functional inorganic materials such as ceramics and metals are used to fill the interconnected pores. The main challenge of these filling techniques discussed earlier is that they lead to the elimination of the porous structure of 3D nanonetworks in a composite state. That is where ALD comes in. The feature of conformal coating of ALD of the surface of porous structures is desirable. Surface properties can be changed, and the pore size can be reduced while we maintain a porous structure. With ALD, it is possible to create a lightweight and highly porous material. It is important to maintain the least possible process temperature below the polymer melting point temperature to preserve the 3D nanostructure during atomic layer deposition of inorganic thin films on 3D nanonetworks [66][67][68]. Figure 3 shows an  i. Gentle deposition process for sensitive substrates ii. Low temperature and stress iii. Coats Teflon iv. Good adhesion • Inherent film quality associated with self-limiting • Self-assembled nature of the ALD mechanism • Depositing multilayers is possible • Stoichiometric control • Relatively low system temperatures illustration of (a) a cycle of thermal/plasma-enhanced atomic layer deposition, (b) the manufacturing process of hollow inorganic materials or polymer-inorganic composite from a polymer substrate with asymmetrical 3D nanonetworks, and (c) the manufacturing process of hollow inorganic materials or polymer-inorganic composite from a polymer substrate with regular 3D nanonetworks.

Atomic layer deposition simulation methodologies
Atomic layer deposition can be studied experimentally or numerically. There are many numerical methodologies employed in atomic layer deposition, among which are Monte Carlo, Knudsen number, molecular dynamics (MD), computational fluid dynamics (CFD), and DFT. In this review, we discuss two of these methodologies: DFT and molecular dynamics.

Density functional theory
DFT is applied in material science to compute the electronic configurations of solids and atoms [69][70][71][72][73]. It is defined as a quantum mechanical (QM) modelling methodology applied in chemistry, material science, and physics to describe the electronic structures of a many-body system (predominantly the ground state), precisely atoms, the condensed states, and molecules [74]. Nowadays, it is used in binding energy calculations of molecules and the band structure of solids in chemistry and physics, respectively [74,75]. In atomic layer deposition, the DFT can be used to explore the reaction mechanism between the first precursor and the substrate. It can be used to calculate the electron density and electronic configurations at this stage. Where the first precursor fails to bond with the surface of the substrate, DFT can be used to explain why this happens. DFT calculations can also be used to improve the conditions of the substrate so that it is favourable for bonding with the precursor. After purging the by-products and excess first precursor, the second precursor is pumped into the reaction chamber. DFT can be used to map the reaction pathway between the first and second precursors. It can be used to predict the nature of the reaction and the by-products from this reaction. DFT can be applied at the end of the ALD process to predict the thin film's mechanical and electrical properties. DFT can be used to validate the experimental conclusions, differentiate between possible opportunities to be explored, or create a foundation for experiments [23]. DFT aims at calculating the electron density and not the wave function of electrons. Functionals are used in this theory to ascertain the many-electron system properties. The functionals in question are of the spatially dependent electron density in DFT [76]. Of paramount importance in DFT calculations and modelling is the Schrödinger equation [77]. The equation can be simply stated as follows: The symbols are defined as follows: Hrepresents the Hamiltonian operator, ψdenotes a collection of Hamiltonian solutions (eigenstates), E -Eigenvalue, according to Sholl and Steckel [78], when several electrons interact with multiple nuclei, the following Schrödinger equation form is used. Figure 3: Illustration of (a) a cycle of thermal/plasma-enhanced atomic layer deposition, (b) the manufacturing process of hollow inorganic materials or polymer-inorganic composite from a polymer substrate with asymmetrical 3D nanonetworks, and (c) the manufacturing process of hollow inorganic materials or polymer-inorganic composite from a polymer substrate with regular 3D nanonetworks [66].
The symbol m represents the electron mass. The first term in bracket defines the kinetic energy of each electron, the second establishes atomic nuclei to electron interaction energy, and the last is electron to electron interaction energy [78]. For the Hamiltonian, ψ defines the electronic wave function. It is a function of each of the N electrons' spatial coordinates. The Hartree product [79,80] can be used to estimate ψ as a product of the individual electron functions. The following is the Hartree equation: The symbol N represents a number of electrons. According to Sholl and Steckel [78], DFT is based on two significant mathematical theorems, which were derived by Kohn and Hohenberg, as well as the Kohn and Sham equations, which were derived in the mid-1960s. Kohn and Hohenberg proved that Schrödinger's equation ground-state energy is a distinctive functional of the electron density [81]. By using the result of Kohn and Hohenberg, we can express the ground state energy E as E[n(r)], where n(r) is the electron density. The major weakness of this theorem is that it does not define a functional. This led to Hohenberg-Kohn's second theorem [82], which addresses an essential property of a functional. It explains that the electron density minimising the energy of the total functional is the true electron density related to the complete solution of the Schrödinger equation [77,78]. This energy functional can be written as follows: Excl i is a function that describes the exchange-correlation functional. The phrase encompasses the quantum mechanical phenomena that are not included in the "known" concepts of quantum mechanics. The known terms come from four contributions as outlined in the following equation: The first right term is the kinetic energies of the electrons, the second is the electrons to nuclei Coulomb interactions, the third is the pairs of electrons Coulomb interactions, and the last is the pairs of nuclei Coulomb interactions. Kohn and Sham proved that finding the right electron density requires a solution of a set of equations in which each equation contains a single electron [82]. The following is the Kohn-Sham equation: where V(r) is a potential that defines the collection of atomic nuclei to electron interaction, V H (r) is a potential describing the Coulomb repulsion between one electron under consideration in one of the Kohn-Sham equations and the total electron density defined by all electrons in the question, and V XC (r) is a potential defining exchange and correlation contribution to the single-electron equations. It is a "functional derivative" of the exchange-correlation energy.
Kohn-Sham equations [81,83] are solved by defining the Hartree potential [79]. Defining the Hartree potential requires us to understand the electron density. It is important to note that finding the electron density requires one to know the wave functions of single electrons, but these can be known when the Kohn-Sham equations are solved. A circle is thus created. Sholl and Steckel [78] outline how somebody may break down this circle by using the following algorithm: 1) A provisional electron density n(r) is defined.
2) The initial provisional electron density is employed in finding the solution (the single-particle wave functions, ψ i (r)) of the Kohn-Sham equations. 3) Next, the electron density, which is defined by the Kohn-Sham single-particle wave functions, is calculated from 2, ( ) ( ) ( ) = ∑ n r ψ r ψ r 2 ks i i i ⁎ . 4) Finally, the computed electron density n KS (r) is compared with the trial electron density n(r), which was used to solve the Kohn-Sham equations. The trial electron density is adopted as the ground-state electron density if the two densities are the same. It can then be used in total energy computation. The provisional electron density must be adjusted if the two densities are different. The process then starts all over again.
A means for the use of the Kohn-Sham equations [82] is offered by the uniform electron gas. The exchange-correlation potential [76,82,84] must be adjusted to be the known exchange-correlation potential from the identical electron gas at the electron density seen at each position. The following local density approximation (LDA) equation is used.
For the Hartree-Fock method [79,80,85], approximating the N electrons wave function, by assuming for a moment that the electrons do not affect each other, in this case, the electron's Hamiltonian is given by: where h i defines the energies (potential and kinetic) of the ith electron.
The Hartree product does not fulfil all the essential criteria for wave functions. The major disadvantage of the Hartree product is that when two electrons are exchanged, and its sign does not change. This draws a need for the use of a superior estimate to the wave function. The Slater determinant is used as an approximation. One-electron wave functions are brought together to form the N-electron wave function, which will satisfy the antisymmetry principle. The antisymmetry principle explains that if two electrons exchange positions, the wave function ought to change the sign because electrons are fermions. The following is the Slater determinant: where 1 2 is a normalisation factor. x i is a coordinate vector defining the ith electron position and its spin state, and χ j (x i ) represents a set of spin orbitals.
For a Hartree-Fock (HF) calculation [79], the atomic nuclei positions are fixed, and then, the N-interacting electrons wave function is determined. Each electron's Schrödinger equation has the following form: where V H (r) is the same Hartree potential. It is defined by the following equation: When solving a single-electron equation practically, the spin orbitals must be defined using a finite amount of information. The finite function sets are described that can be summed up for approximating the exact spin orbitals. The spin orbitals are approximated as follows: where α i,j is the expansion coefficient and ϕ i (x) is the basic set for the computation (where i = 1,2,3,…,k).
Defining the Hartree potential in the equations of single electrons requires that the electron density is known. Knowing the electron density needs prior defining of the electron wave function (which is computed using the separate spin orbitals) [78]. A Hartree-Fock computation is an iterative procedure conducted as follows: through the specification of expansion coefficients, α j,i . 2) Define the electron density n(r′) from the approximate spin orbitals. 3) Solve the spin orbitals equations using the answer from step 2. 4) The solutions are adopted as the solutions to the Hartree-Fock problem if the spin orbitals found in the third step are the same as those used in the second step. If not consistent, a new estimate should be used, and then, the procedure is restarted from the second step.
For the best results, the DFT methodology is applied at the atomic scale. DFT has proved to be useful in the forecasting of molecular properties as well [86]. It can be used as a counterpart of experimental studies. In the ideal situation, the calculated characteristics are accurate enough to distinguish between structural replacements for chemical intermediates or other species that are not sensitive to the explication of experimental structure cdd [86].

Applications of DFT calculations in atomic layer deposition
When ground state DFT calculations are performed, one can get various routine quantities such as equilibrium geometry at absolute zero, equilibrium lattice structure, formation energy, partial charge at each atom, magnetic state of the system, fermi level, band gap, optical absorption spectrum, and spin properties [87][88][89][90][91]. In equilibrium geometry, the calculation yields the geometry with the minimum forces on each of the atoms in the system. Band gap calculations tell us the nature of the material, for example, if it is metallic, semiconductor, or an insulator. Metallic materials have a band gap of 0 eV, while insulators have a very large band gap of approximately 15 eV.
In surface science and catalysis, DFT calculations can be useful in obtaining the surface properties of a material. These include surface energy, the adsorption energy of a gas molecule on a surface, adsorption configuration of a gas molecule on a surface, bond distance, charge transfer, reaction energy barrier, and work function. DFT is useful in predicting the adsorption energies and configuration of how the gas bonds with the surface, for example, if the gas is CO, DFT calculations can predict which atom between oxygen and carbon bonds with the surface. In charge transfer, DFT can predict which of the two materials is the charge donor and which one is the charge acceptor among the molecule and the surface.
In semiconductors, DFT can predict the electronic band structure, defect bands, conduction bands, valence bands, and band alignment across interfaces [92][93][94]. The band structure is used to explain the electrical properties of the material. DFT can predict the nature, location, and type of the defect.
DFT application is superconductors is still in its infancy. DFT calculations can predict electron-phonon coupling and superconducting critical temperature, which is sometimes referred to as the critical transition temperature. The application of DFT to predict electron-phonon coupling has yielded great success, but the same cannot be said for predicting critical superconducting temperature.
In molecular electronics, DFT can be used to predict conductance as a function of the energy and current-voltage characteristics of a material.
In this article, we discuss the applications of DFT in surface science and catalysis and also in semiconductors.

Applications in surface science and catalysis problems
DFT can be used to investigate area-selective atomic layer deposition (ASALD) of SiO 2 in an ABCtype cycle with acetylacetone used as a chemoselective inhibitor [52]. In this cycle, step A entails exposing the surface of the substrate to inhibitor molecules, whereas steps B and C are the usual conventional ALD process steps. Inhibitor molecules deactivate certain areas of the substrate and prevent the growth of films on those regions during the ALD cycles (see Figure 4). ASALD has a great potential to be applied in self-aligned manufacturing systems for future-generation nanoelectronics [52,95]. The function, projector augmented wave (PAW) [97,98] was used by Mameli et al. [52] as executed in Vienna ab initio simulation package (VASP) to perform electronic structure calculations [98][99][100]. They carried out first-principle computations using the generalised gradient approximation (GGA) to the DFT. Two distinguishing binding orientations were identified by their DFT calculations for physisorption.
Gas-phase density functional theory (GPDFT calculations can be used to explain valence and core levels gasphase photoelectron spectroscopy of transition metal alkyla-mino complexes, M(N(CH 3 ) 2 ) 4 . where M in this case, is Titanium (Ti) or Hafnium (Hf) [38]. These calculations can be used to predict atomic layer deposition reactions on the surfaces of these molecules. In their research, Shayesteh et al. performed DFT calculations using the hybrid M06-2X [101] functional as executed in the Gaussian 09 code [102]. The double-ζ 6-31G(d,p) basis set was used for the hydrogen, carbon, and nitrogen atoms. The LanL2DZ basis set [103][104][105], together with the Los Alamos effective core potential, was utilised for the hafnium and titanium atoms. Population analysis techniques of Kohn-Sham orbitals were used to allocate valence spectra of Hf(N(CH 3 ) 2 ) 4 and Ti(N(CH 3 ) 2 ) 4 compounds. Through the use of DFT calculations, computed vertical ionisation energies and adiabatic ionisation energies were matched with their data from experiments.
DFT, coupled with Monte Carlo simulations, can be employed to model equilibrium conformation of Pt 2 Ru 3 nanoparticles [106] utilised in polymer electrolyte fuel cells as anode catalysts (PEFCs) [107][108][109]. Alam, Saito, and Takaba conducted DFT calculations to approximate (using slab configurations) the binding energy. Pt-Ru alloys significantly enhance the catalytic performance in the electrochemical oxidation of carbon monoxide from carbon monoxide-contaminated H 2 fuels or suppression of carbon monoxide adsorption on the catalyst [110][111][112].
The DFT calculations showed that in the Pt 2 Ru 3 alloy slab configuration, platinum atoms covered the top surface, or a more significant percentage of platinum led to higher stability compared to other conformities of atoms. The second layer is favourably comprised ruthenium atoms or a greater proportion of ruthenium atoms for stable structures. The quantity of Pt-Ru bonds in Pt 2 Ru 3 increases as the particle size also increases. Figure 6 shows the different slab configurations.
DFT was applied in ALD-triggered iron-indium-sulphur cluster and gradient energy band in zinc-indiumsulphur photoanode to improve the reaction of oxygen evolution [113]. Meng et al. did their DFT calculations using VASP [113]. The team designated Perdew-Burke-Ernzerhof (PBE) functional for the electron interactions of the exchange correlation [114]. Meng et al. [113] chose the projected augmented method [97] to hypothesise the pseudopotential for treating the interactions of ion to electron [97]. For all other calculations, the team employed panewave cut-off energy of 400 eV and 4 × 4 × 1 Γ-centred Monhorst-Pack k-point mesh [115]. DFT-D3 method of Grimme was adopted to define the interactions of van der Waals [116]. DFT has also been used to investigate exceptionally dry methane reforming catalysts with enhanced in situ produced Ni-Fe nanoparticles on perovskite through ALD, which have been demonstrated to be extremely effective [117]. Similarly, some studies aimed to develop active and durable catalysts [51,[118][119][120]. DFT calculations were used to reveal the method of the improved exsolution via atomic layer deposition. The DFT calculation results showed that the order of Ti < Cr < Fe < Mn < Co < Ni < Cu was followed by the co-segregation energy. This implied that bulk Ni could be replaced with deposited iron (Fe) [117]. Figure 7 shows DFT-computed cation exchange energies of deposited bulk Ni and transition metals. The deposited metal can freely exchange with Ni as indicated by the negative sign in exchange energy.
Density functional theory has been utilised to investigate the growth mechanism of zinc oxide ALD [121]. In their investigation, Afshar and Cadien [121] used GAUSSIAN 09 [102] to execute molecular orbital calculations. To determine the geometry of the stationary points, the corrected density functional approach (B3LYP gradient) was used [70,121]. In their conclusion, Afshar and Cadien [121] suggested a comprehensive atomistic mechanism for zinc oxide ALD growth employing the DFT study of (C 2 H 5 ) 2 Z n and H 2 O half-reactions on a zinc oxide surface. By using this DFT approach, they showed that the two reactions proceeded through the creation of complexes on the surface. The first growth mechanism of zinc oxide ALD on hydroxylated Si was investigated by Ren [122]. He presented the impact of single and double hydroxyl sites on the initial phases of the atomic layer deposition growth [123].
DFT was used in the preliminary reaction mechanism investigation of atomic layer deposited Pt on a hydroxylated graphene surface [19]. In their investigation, Mohlala et al. used molecular oxygen and methylcyclopentadienyl trimethylplatinum (MeCpPtMe 3 ) as their precursors. In their research, the reaction pathways were analysed using DFT calculations under Dmol3 code in Materials Studio software, a product of Accelrys Software Inc [124]. Their DFT calculations established that the elimination reaction of the methyl group bonded to platinum happens more easily than the freeing of the MeCp ligand [19]. They adopted two possible configurations to investigate which interaction was more stable between the precursor and the site of adsorption. Models (a) and (b) were created in which model (a) interacted through the MeCp ring and model (b) interacted through one of the methyl groups, which is bonded to platinum. Their DFT simulation and calculations showed that model (b) is the most stable configuration [19]. Figure 8 shows the preliminary orientation of MeCpPtMe 3 on the surface of grapheme, which was hydroxylated.
Graphene is now being used in the development of next-generation functional nanoparticles for numerous applications. DFT was employed in finding an estimated Schrödinger equation solution. They performed a DFT simulation of infrared (IR) characteristics in support of the proposed mechanism through a broad clarification of forming and breaking of bonds during each phase of the reaction. Figure 9 shows the structures, transition states, intermediate states, and dissociation reaction products. Figure 10 shows the infrared data when graphene oxide matrix adsorbed Zn-H, which was computed using DFT, and the associated optimised structures for several Zn 2+ adsorption steps [77].
Mohlala et al. [128] employed DFT in their article on the thermal stability of Ti halide precursors for the TiO 2 ALD on Pt(111) substrate. DFT was employed to ascertain the reactivity and thermal stability of Ti-based halide precursors. They also used DFT in conjunction with the computerised numerical technique and the dissociation of bond ligand automation tool from the Schrödinger material science package [129][130][131] to study and evaluate the thermal stability and reactivity of the precursors. The reason for the use of DFT is that they knew it as a promising technique to understand the chemistry of innovative ALD precursors and enabling precursor development, which has high selectivity and reactivity and thermal stability [132][133][134]. By using DFT simulation and calculation, Mohlala et al. predicted Ti tetrahalide (fluorine, chlorine, and iodine) precursors through the use of Schrödinger's automated enumeration tool. They predicted that a library of three halides forms 81 different complexes, 15 of which have exceptional ligand compositions [128]. The precursors were improved utilising DFT as realised in Jaguar [131]. Figure 11 shows an illustration of 15 of the precursors.
Becker and Sierkar [70] employed DFT in their atomistic simulations research of plasma-enhanced ALD. Their primary investigation of the precursor products with hydroxylated silicon dioxide surfaces used DFT calculations and showed that reaction energies were highly dependent on the surface structure [20], specifically for the precursor which is bonded to two oxygen atoms on the surface.  Time-dependent DFT in conjunction with the Ehrenfest dynamics scheme was employed by Lee et al. [135] for the investigation of the dissociation of cobalt tricarbonyl nitrosyl (Co(CO) 3 NO) ligands caused by electronic excitation. The team scrutinised excited-state non-adiabatic dynamics of ligand dissociation of Co(CO) 3 NO. Cobalt tricarbonyl nitrosyl is a precursor mainly utilised to grow cobalt thin films in innovative technologies, where the reaction of the precursor is improved by electronic excitation. They studied the Co (CO) 3 NO microscopic details of dissociation motivated via electronic excitation and established two direct and distinct dissociation paths of the NO ligand on the cobalt tricarbonyl nitrosyl (Co(CO) 3 NO).
Murray and Elliott [134] employed DFT in the prediction of the composition of ALD-grown ternary oxides. They calculated using DFT methods the geometry optimised structures and energies for 17 metal precursors with five different anionic ligands and their corresponding hydrolysed complexes. The hydrolysis DFT models were utilised to effectively clarify the reaction and stoichiometric result with regards to metal cation ratios, which were seen through experiments for different ALD-grown systems for ternary oxides [134]. When it comes to electronic interchange and correlation, transition metal complex computations are prone to inaccuracies; therefore, the selection of the proper DFT functional is vital to yield correct geometries and enthalpies of reactions. Their DFT calculations confirmed that metal complexes that contain the alkoxide O t Bu and β-diketonate thd ligand are substantially unreactive with regard to the removal of H-thd during atomic layer deposition. Most of the ALD ternary oxides were successfully predicted by the DFT models that were used in their research for the interactions of metal precursors with oxide surfaces, demonstrating that DFT models may be used to predict experimental metal cation concentrations.
Dkhissi et al. used DFT in multiscale modelling of HfO 2 ALD of thin films, which were grown on silicon [136]. When the DFT model was applied in the hybrid functional, it was possible to determine the complete physicochemical processes and their related energies of the two half-cycles that occurred at the initial stage of the film development. It was demonstrated in their work that a novel kinetic Monte Carlo (KMC) methodology could be developed within a multiscale approach, starting with molecule-surface interactions that were primarily treated at the DFT level of  modelling and progressing to atomic-scale film growth that was conducted using a KMC technique [136].

Applications in semiconductors
Weckman et al. [91] utilised DFT in the study of zinc oxide atomic layer deposition. In their research, they investigated the ZnO thin film growth, which was grown via DEZ/H 2 O process. To achieve this, they modelled the surface chemistry using first-principles kinetic Monte Carlo (KMC). The KMC enabled them to easily extract data equivalent to experimental results by incorporating DFT calculations done on the ZnO(100) surface into a kinetic model. ZnO is a semiconductor, and its properties include high transparency, piezoelectrical properties, and tunable electrical conductivity. It can be applied in transistors, sensors, and solar cells [137][138][139]. The results of the DFT-simulated mass increment data during an atomic layer deposition cycle were found to tally with the experimental findings and showed that the majority of the deposited layer is formed during the DEZ pulse [91]. The activation energy for chemical processes was also calculated using DFT.
Sun et al. [140] used first principles DFT calculations to study the chlorine residue effects on the ALD of Hafnium oxide (HfO 2 ). HfO 2 -based gate dielectrics are thermodynamically stable with silicon. The density of the state was calculated using the CASTEP code utilising the local density approximation method to get the electronic band structure of HfO 2 . Their calculations yielded a band gap of 6.69 eV, which was comparable to the reported experimental value of 5.6-5.9 eV since they had neglected excitonic and phonon broadening effects, which typically would have reduced their calculated value by 0.8 eV. The team noted that chlorine decreases the band gap of hafnium oxide (by 0.58 eV) by increasing the valence band maximum (VBM). There are two types of defects in hafnium oxide (substitutional and interstitial defects), and these two harm the hafnium oxide gate dielectric electrical performance [140].
Kelaidis et al. [94] used DFT to study the Sn 1−x Pb x O electronic properties and the SnO/PbO system band alignment. The SnO grown by atomic layer deposition was doped with Pb to form Sn 1−x Pb x O. The team conducted DFT calculations to investigate the electronic and structural properties of Sn 1−x Pb x O. They investigated the orbital effects in the valence and conduction band states. They also investigated the effect of the intrinsic defects on band alignment on Sn vacancies in SnO and Oxygen vacancies in PbO. The team found out that Pb doping provided a way of SnO band gap tuning without the introduction of defect states in the gap region. They also discovered that a 0.4 eV increase of SnO was observed when 2% Sn vacancies are introduced. The team noted that the energy band gap of Sn 1−x Pb x O increased virtually linearly from 0.4 up to 1.4 eV as the x value is increased from 0 to 1. Figure 12 shows the computed band gap with regards to the Pb content in the compound.

DFT calculation challenges
The Schrödinger equation cannot be solved precisely for a multi-body system; hence, DFT simulations use approximations. The accuracy of the adopted approximation governs the usefulness of their DFT simulations for several properties [141]. The precise handling of electron correlation founds a significant challenge in DFT. Electron correlation is frequently described with regard to a Hartree-Fock wave function. According to the variational principle, the lowest energy orbitals may be found by leveraging the HF wave function, a single-configuration wave function derived by employing the variational principle. The difference between a system's exact energy and its HF energy at the entire basis set limit is typically used to calculate electron correlation energy [141].
The most frequently utilised DFT form is the Kohn-Sham formulation [142]. The KS-DFT signifies the density using a single Slater determinant and has a mean-field formulation, but it is an exact theory [141]. The energy in Kohn-Sham density functional theory is not computed by additional configuration state functions (CSFs), but by employing an exchange-correlation functional, that is a density Figure 12: Computed band gap with regards to the Pb content in the compound [94]. functional. We must, therefore, estimate the functional of the exchange correlation. The existing functionals are more precise for correlated systems that are weak than the strong ones. Self-consistent-field variational computations are used to determine the orbitals of a Kohn-Sham density functional theory determinant. The variational principle at times may account for static correlation by using a slater determinant with different symmetry properties from the precise state wave function.
KS-DFT can lead to an error source of delocalisation, sometimes called self-interaction error [143]. By substituting nonlocal high-frequency exchange for a portion of the local exchange-correlation functional, the self-interaction error can be reduced [85,144]. DFT does not treat discrete electrons in precisely identical ways but uses their total density. This means any one-electron system does not play a distinct unique role in DFT. The Coulomb energy for one electron is cancelled by the exchange energy [145]. The orbital rotation operator that is used for generating a foundation for time-dependent DFT equations generates only the configurations that are singly excited with regard to the ground state. This means higher order excitation states (e.g., the excitations with a high degree of double excitation character) are excluded in time-dependent DFT. However, they can be incorporated through the exchange-correlation functional [142].
Many variants of the modern density-built subsystem and embedding approaches [146][147][148] use potential reconstruction techniques [149][150][151]. This leads to a target density being rationalised numerous times in a self-consistent loop. The potential matching that target density needs to be calculated at every step of the loop, leading to a considerable added computational expense to simulations that make no use of effective inversion methods.
The other challenge of DFT is its complexity. DFT loses its simplicity when DFT functionals become complicated like the full configuration interaction (FCI). This is factual of its nature of computation. However, the attractiveness and challenges of DFT lie in between simple and very complex [145].

Molecular dynamics
Molecular dynamics (MD) simulation [70,[152][153][154] is a computational procedure for computing the transport properties and equilibrium of a classical many-body system [155]. MD focuses on the dynamics of atoms, molecules, and clusters in the condensed and liquid states [156]. It is used in many engineering and science disciplines like chemistry, physics, and materials engineering to determine the motion and equilibrium conditions of each atom or molecule.
During the first stage of the atomic layer deposition process where the first precursor is forced into the reaction chamber, molecular dynamics calculations can be used to predict the kinetic properties of molecules such as material structure and complex and the stability of the substrateprecursor 1 bond. It can predict the binding energy between the substrate and the first precursor molecules. In the second stage where the by-products and excess first precursor are purged, MD can be applied to predict the mechanism of removal. In the third stage where the second precursor is pumped into the reaction chamber, MD simulations can be used to predict the nature of the reaction, the binding energy between the first and second precursor atoms, and the stability of the bond. It can also be applied in the fourth stage the same way as the third stage since these purging stages have the same principles. Molecular dynamics can be used to predict the thermodynamic properties of the deposited thin-film material, based on intermolecular interactions. The deposited thin film can be analysed by MD simulations to predict free binding energy, the heat of vapourisation, and boiling point. According to Sator [157], the molecular dynamics method is valid in both equilibrium and non-equilibrium physical phenomena. Molecular dynamics is a useful computational approach that may be used in simulating various physical phenomena [158,159] if there is satisfactory computation power. Besides ReaxFF [70,[160][161][162], other computational software that can be used in molecular dynamics simulations are LAMMPS [163][164][165][166], NAMD [167], GROMACS, CHARMM [167], and AMBER [168,169]. The molecular dynamics method uses Newton's equations of motion to describe the motion of molecules. It is simple and logical [170]. Simulation of particle flow on a computer MD program follows the equations of motion [157,171]. The force exerted by an external field and ambient molecules is defined by the following equation of motion: where the mass of the molecule I is denoted by m i .
The Taylor series expansion can be used to transform equation (13) into an algebraic equation through expanding the second term differential term into a Taylor algebraic expression as follows: Another Taylor series expansion form is required to estimate the second-order differential term in equation (13). This is:

(15)
The second-order differential term may be solved by removing the first-order differential term from equations (14) and (15) as follows: The right-hand side last term indicates the approximation accuracy. Terms higher than h 2 are neglected in this case. The expression is called the central difference approximation if equation (16) is approximated as follows: The x-component of Newton's equation of motion may be stated by using this approximation with the notation r i = (x i , y i , z i ) for the molecular location and f i = (f xi , f yi , f zi ) for the force acting on particle I as follows [157]: Equation (18) is an algebraic equation; therefore, the next time step position of molecules can be calculated using the current position, past positions, and the current force. Equation (18) does not need velocity to determine the next time step position of molecules. If velocity is required, it is calculated from the central difference approximation. The velocity equation is expressed as follows: Noting that the positional first and second-order differentials are identical to the velocity and acceleration, in that order, ignoring differential terms equivalent to or higher than the third order in equation (14) yields: Since the term of velocity appears on the right, another equation to specify velocity can be employed.
The force term can be modified to improve accuracy. The equation becomes: When the central difference approximation is applied to the first-order differentials, we get these two equations: According to Sator [157], when using the velocity Verlet to carry out a molecular dynamics simulation, the primary procedure is as follows: 1. The initial position and velocity of the molecules must be hypothesised. 2. Next, the forces that act on molecules must be calculated. 3. The next time step positions of molecules must be evaluated from equation (20). 4. The next time step velocities of molecules must be evaluated from equation (22). 5. Procedures from the second step are repeated.

Applications of molecular dynamics simulations
Molecular dynamics simulations can be used to predict the kinetic properties of molecules [95,158,169,172]. Kinetic properties are time-based evolutions of conformations. These include refinement of material structure and complex, stability of structures, conformational sampling of small molecules, protein folding, stability of protein-ligand interactions, and transport of ions and ligands.
MD can be used to predict the thermodynamic properties of a material. These predictions are based on intermolecular interactions. MD simulations can predict free binding energy, heat of vapourisation, pressure, boiling point, Poisson-Boltzmann or generalised Born and surface area continuum solvation, and free energy perturbation (FEP).
Molecular dynamics is most suitable for the application at the nanoscale. There are various applications of molecular dynamics as identified by Gunsteren and Berendsen [173]. The following are some of the applications: i. Interpretation of biochemical [174,175] data on repressor-operator binding. ii. Interpretation of biophysical [176] data on membrane properties. iii. Refinement of structures based on nuclear magnetic resonance (NMR) data [112,177]. Tampieri defines NMR spectroscopy as an analytical method that provides data on the local magnetic field around atomic nuclei [178]. Its principles are that the structural or chemical composition of substances can be known from their nuclei, which have their distinctive magnetic field. iv. Refinement of structure based on crystallographic Xray or deutron diffraction data [179]. X-ray crystallography is used to determine three-dimensional crystalline structures at the atomic scale [180]. v. Prediction of structural changes. vi. Prediction of free energy changes.

Applications of MD in predicting the kinetic properties of molecules
Molecular dynamics simulations can be utilised in the atomistic simulations study of plasma-enhanced ALD [20,47,70]. In their research article, Becker and Sierka [70] presented a simulation structure for plasma-enhanced atomic layer deposition processes by compounding the Monte Carlo [157,179] deposition algorithm and structure relaxation using MD to account for the steric interruption and the surface overlap restrictions consistent with the actual step of precursor deposition. The structure was relaxed after each step of attachment, using either local structure optimisation or short-simulated annealing using MD. They calculated the chemisorption energy and used it to compute the acceptance criterion using the Metropolis-Hastings rule [181]. They also applied their simulation procedure to study the PEALD synthesis of silicon dioxide thin films employing a bis-diethylaminosilane precursor. Preliminary examination of the precursor products with hydroxylated silicon dioxide surfaces using DFT calculations showed a strong reliance of reaction energies on the structure of the surface, specifically for the precursor that was bonded to two oxygen atoms at the surface. Their PEALD simulations showed that the precursor binding to one oxygen atom at the surface favours amorphous layer growth, many -OH impurities, and the voids formation. Molecular dynamics simulation can be applied to investigate the sinking of graphene during chemical vapour deposition (CVD) growth on a substrate of semi-molten copper [119]. The principles of the molecular layer deposition procedure are synonymous with ALD since, in both processes, the films are created by successive, self-restraining surface reactions [182][183][184]. Xu et al. executed systematic MD simulations to discover graphene nanostructures stability (graphene nanoribbons and carbon nanoclusters) on a substrate of semi-molten copper. Figure 13 shows C 24 sinking on the surface of Cu(111) at various temperatures.
Molecular simulations and Markov state modelling can be used to disclose the structural variety and dynamics of a theophylline-binding ribonucleic acid aptamer in its unbound state [185]. Warfield and Anderson applied MD simulations to characterise the ensemble of conformations adopted by the aptamer in the absence of theophylline and also applied the Markov state modelling to forecast the transition kinetics between liberated conformational states. Ribonucleic acid aptamers are oligonucleotides that connect with high specificity and attraction to targeted ligands [185]. It was observed in their simulations, as well as in the experiments, that the theophylline binding site may be located in a variety of binding-incompetent states, all of which lack a binding pocket capable of holding theophylline. Binding-incompetent states can become binding competent by rearranging the binding site structurally on a period ranging from nanoseconds to microseconds.
MD can be used to simulate the film structure of Al 2 O 3 and its composition during atomic layer deposition [186,187]. MD simulations allowed them to trace the progression of the Al 2 O 3 film microscopic structure as it grew in thickness [186]. They found that the surface of Al 2 O 3 is ended mainly by hydrogen and oxygen since they are the favourite for the ambient. This enabled them to conclude why hydrogen impurities in atomic layer-deposited Al 2 O 3 films are generally low.
Zheng et al. [10] used MD simulations in their research on modelling and in situ probing of surface reactions in ALD. As part of their investigation into the effect of early surface chemical states on the reaction kinetics and the creation of traps during ALD nucleation, they looked at Al 2 O 3 ALD deposition kinetics [148] on hydrogenated and oxidised Ge surfaces. Their simulations included ReaxFF [70,154,160,161], molecular dynamics (MD), and nudged elastic band (NEB) [188] to mimic the deposition process, including adsorption and breakdown of ALD precursors, as well as the diffusion and desorption of their products produced. Their simulations accurately predicted the nucleation features, which were in close accord with the available literature and experiments. Figure 15 shows the TMA and H 2 O dose in their ReaxFF simulation.
Xiao et al. [154] used MD simulations to investigate the surface chemistry of Cu metal and CuO ALD from copper(II) acetylacetonate. The results of this simulation showed copper(II) acetylacetonate chemisorbs on the hollow site of the copper(110) surface and decays readily into a copper atom and the acetylacetonate ligands. Figure 16 shows snapshots for the reactive molecular dynamics dissociation of copper(II) acetylacetonate on the Cu(110) surface. Figure 17 shows the reactive MD snapshots for the reaction between copper(II) acetylacetonate and hydrogen atoms on the copper (110) surface. Atomic hydrogen was highly reactive towards the copper precursor. The Cu (acac) 2 broke the Cu−O bonds upon hydrogen collision forming a H 2 (acac) molecule. It was then released to the  [154].
Ab initio MD simulation of a copper monolayer aggregation on a WN (001) surface was investigated by Bo et al. [189]. They carried out an ab initio molecular dynamics simulation with a copper monolayer at 500 K on the surface, which was completely equilibrated. Figure 18 shows the molecular dynamics chosen snapshots of the trajectory. The copper monolayer was initially well aligned and in proportion with the substrate. Through surface diffusion, copper atoms were then rapidly disassembled and moved away from their equilibrium locations when molecular dynamics was run. At t = 1 ps, some copper atoms began to move upward and grow three dimensionally. At t = 3 ps, the copper atoms were clustered together and formed islands at the surface. An ordered monolayer structure became remarkably disordered. Further molecular dynamics runs prompted copper atoms to constantly move, forming surface islands with different shapes. No wetting of copper clusters on the surface was detected.
Guoyong et al. [190] applied MD in surface pseudorotation in Lewis-base-catalysed SiO 2 ALD. The pathways of the reaction for the uncatalysed and catalysed SiO 2 atomic layer deposition were explored using static transition-state searches and Born-Oppenheimer MD simulations. The team discovered that the reaction for uncatalysed SiO 2 atomic layer deposition undergoes a step-by-step pathway. The step that determines the rate is the Si-O bond development complemented by the rotation of SiCl 4 and has a comparatively high activation free energy barrier.

Applications of MD in predicting a material's thermodynamic properties
Guo and Wang, in their research of thin films of noble metal alloys by ALD and rapid Joule heating [191], performed MD simulations using large-scale atomistic/molecular massively parallel simulator (LAMMPS) open-source code [192,193]. They adopted the embedded atom method (EAM) potentials [194] in portraying the interatomic interactions, which define the lattice parameter, surface energy, cohesive energy, elastic constant, transitional and thermodynamic properties of metals, as well as their alloys. Figure 19 illustrates a schematic of the fabrication process of the thin film alloy of Rh/Ir. From their research and simulations, Guo and Wang concluded that the ALD technique can manufacture fairly uniform and continuous sub-fifty-nanometre thick alloy thin film on 2D substrates when compared to other methods of manufacturing alloys. They also employed MD simulations to explore the melting of the Rh/Ir thin film alloy.
Molecular dynamics can be applied to investigate temperature and nitrogen-to-aluminium ratio effects in the heteroepitaxial growth of aluminium nitride. Zhang et al. [195] noted in their simulations that when the temperature increased from 1,000 to 2,000 K with a nitrogento-aluminium flux ratio of 2.0, the growth rate of the Figure 14: (a) Blue line shows the total density profiles at 0 ALD cycles, red for four ALD cycles, and green for eight ALD cycles with respect to the Z-direction distance of the Al 2 O 3 film at 300°C [186]. (b) The number density profiles for oxygen (blue), aluminium (green), and hydrogen (red) as a function of the Z-direction distance of Al 2 O 3 film at 300°C [186].
aluminium nitride film decreased. They also observed that deposited aluminium nitride crystallinity was noticeably enhanced as the temperature increased between the 1,000 and 1,800 K temperature range, and it became saturated between 1,800 and 2,000 K. Stoichiometry was found to be closely associated with the deposited films crystallinity. Good crystallinity films are related to a near 50% nitrogen fraction. In their simulations, the aluminium nitride growth on the c-plane aluminium nitride was studied by the LAMMPS [193] with the Tersoff potential [196].  shows the fraction of nitrogen of deposited films at a temperature of 1,800 K with various nitrogen-to-aluminium flux ratios starting at 0.8-2.8.
Molecular dynamics simulations can be utilised to investigate the melting of iron nanoparticles with or without defects. Karpina et al. [137] used MD simulations to investigate the thermal properties of bulk iron material and iron nanoparticles (FNP) by means of a ReaxFF reactive force field. The melting behaviours of FNPs from 300 to 2,500 K were chosen to study thermodynamic and energy properties, which include potential energy plots, radial distribution function, and Lindemann index [120]. They observed that the ReaxFF force field could detect the size effect in FNP melting regardless of energy or structure evolution aspect. The following table shows the FNP model information for melting simulations. In their article, Sun et al. state that there are four methods to evaluate materials melting points [197][198][199] in molecular dynamics. (1) Plots of potential energy; (2) system heat capacity C v , and (3) Lindemann index [200] that vary with temperature and equilibrium temperature [161] of the solid-liquid interface (see Table 4). Figure 21 shows one of the ways of evaluating melting points, the C v (heat capacities) of various size FNPs as a function of temperature obtained from simulations.

Molecular dynamics simulation challenges
Just like any other modelling and simulation methodologies, molecular dynamics has its shortcomings. This is discussed in the following paragraphs.    Figure 18: MD trajectory [189]. The copper monolayer was initially well aligned and in proportion with the substrate. Through surface diffusion, copper atoms were then rapidly disassembled and moved away from their equilibrium locations when molecular dynamics was run. Figure 19: Schematic of the fabrication process of the thin film alloy of Rh/Ir [191].
MD simulations can be applied to analyse chemical processes, which include biomolecules allosteric transitions, phase transformation nucleation events, chemical reactions, solids fracture and crack propagation, crystals defect dynamics, and many other applications. These simulations, however, may take several days depending on the size of molecules under analysis [201]. The atom (or molecules) quantity in the system determines the time scale available to the simulation problem, hence the time-scale problem severity [174,202]. This governs the speed at which investigations can be conducted. The time-scale problem is the major reason why brute-force molecular dynamics simulations are not practical when investigating crystal nucleation.
Virtually every computer simulation of the nucleation of crystals in liquids uses some classical nucleation theory (CNT) [202,203] elements. Classical nucleation theory was invented to define supersaturated vapours condensation into the liquid phase although the majority of the concepts can likewise be employed in the crystallisation of supercooled liquids and supersaturated solutions. It is assumed that classical molecular dynamics simulations can handle up to approximately 10 5 molecules on a nanoseconds or microseconds time scale. There is a very thin set of situations for which brute-force classical molecular dynamics simulations may be applied to explore nucleation, and it is only typically at strong supercooling.
Most MD simulation packages are robust and need adequate computing power [204]. They also require a substantial amount of time for one to familiarise themselves with them. Experiments are still far much superior to simulations [202]. There have been instances where molecular dynamics simulations have not adequately predicted the outcome of experiments. This makes it necessary to couple simulations with experiments.
Molecular dynamics simulations generate large amounts of data, which intrinsically presents big data challenges [205]. Storage, dissemination, and management of these data in the range of terabytes are significant even though the capacity of storage resources is increasing [206].
It is difficult to fully capitalise on the richness of data that is offered by simulations where each molecule is indicated in atomistic detail. Simulations are usually done with a definite scientific problem in mind, and we often concentrate on that problem in mind during the investigation in a stringently hypothesis-driven fashion. Data that seem not to answer the question in mind is usually neglected. This hinders discoveries that could come from the neglected data.
For molecular dynamics simulation results and experiments to agree, they must be carried out under the same conditions. This may not always be feasible since simulations are usually carried out under ideal conditions. Several variables may be present when conducting experiments.
The reliability of the simulation software used can pose a challenge to amateurs and people who have recently started the molecular dynamics field of study [35]. One needs to be thoroughly vexed with all the available software, their application, strengths, and weaknesses. The simulation software quality relies on the attention with which it was built and tested by its developers and their understanding of the molecular models and algorithms incorporated in their software. Testing of simulation packages is usually done on numerous levels within diverse communities of researchers [35].

DFT and MD use by non-specialists
DFT and molecular dynamics simulations use complex software that require adequate time for one to learn. A general appreciation of DFT and MD methodologies is not adequate for one to start using these simulation tools, as this may lead to wrong results. Users of these methodologies must have a coding background, or must at least learn a few coding languages such as python. There are various tutorials accessible from the Internet where one can get started. Non-specialists should seek assistance from DFT and MD specialists if they have to apply these tools in their work.

Conclusion
Computational modelling and simulation methodologies are very powerful tools for analysing microstructures and nanostructures. They complement ALD experiments but have the advantage of significantly reducing the financial burden of doing experiments. Their other benefit is that they can be used to analyse very complex structures and reactions that may be otherwise very dangerous to conduct experimentally in the lab. Through simulations, any problem can be interrogated regardless of practical limitations inhibiting experiments, including hypothetical situations that could not be realised through experiments under any conditions. However, it should be noted that computer simulations remain a theoretical approach whose value lies when new hypotheses and predictions generated can then be subjected to experimental analysis and validation. DFT can be applied in molecule binding energy calculations in chemistry and the solids band structure in physics, for computing the electron density, and for determining the many-electron system properties. DFT can answer problems of understanding material properties in surface science and catalysis, semiconductors, superconductors, and molecular electronics. Molecular dynamics can be used to interpret biochemical data on repressoroperator binding, interpret biophysical data on membrane properties, refine structures based on nuclear magnetic resonance, refine structures based on crystallographic X-ray or neutron diffraction data, and predict structural and free energy changes. Molecular dynamics can predict the kinetic and thermodynamic properties of materials. Several papers have been published where these simulation techniques managed to predict the outcome of experiments successfully. One must be able to distinguish between the discussed methodologies and pick the relevant one for their application. These methodologies can be used to prove and validate results from experiments and predict the outcome from experiments. Most MD and DFT simulation packages are robust and need adequate computing power. Substantial time is required for one to learn the software and produce accurate simulation results. Since the simulation results do not always agree with experimental results, it is important to couple simulations with experiments where possible. Computational modelling and simulations successes greatly outweigh their challenges.
Funding information: The authors would like to appreciate the funding from the Global Excellence and Stature (GES) and University Research Committee (URC) of University of Johannesburg.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
The authors state no conflict of interest.