Investigations on Forming Ether Coated Iron Nanoparticle Materials by First-Principle Calculations and Molecular Dynamic Simulations

The mechanism of coating effects between ether molecules and iron (Fe) nanoparticles was generally estimated using first-principle calculations and molecular dynamic (MD) simulations coupling with Fe (110) crystal layers and sphere models. In the present work, the optimized adsorption site and its energy were confirmed. The single sphere model in MD simulations was studied for typical adsorption behaviors, and the double sphere model was built to be more focused on the gap impact between two particles. In those obtained results, it is demonstrated that ether molecules were prone to be adsorbed on the long bridge site of the Fe (110) crystal while comparing with other potential sites. Although the coating was not completely uniform at early stages, the formation of ether layer ended up being equilibrated finally. Accompanied with charge transfer, those coated ether molecules exerted much binding force on the shell Fe atoms. Additionally, when free ether molecules were close to the gap between two nanoparticles, they were found to come under double adsorption effects. Although this effect might not be sufficient to keep them adsorbed, the movement of these ether molecules were hindered to some extent.


Introduction
Nanomaterials refer to those material which have at least one dimension of structure in nanoscale. They are believed as one the most potential research field all over the world [1]. Taking environmental application as an example, many development opportunities of novel functional materials, processes, and devices are offered by nanomaterials and nanotechnologies [2]. Typically, Fe nanoparticles are well performed nanoscale environmental materials that have been applied for a long time [3]. Besides, Fe nanoparticles are also appropriate for magnetic, electrical, catalytic and biomedical applications [4]. Until recently, the latest research of Fe nanoparticles is still continuing. Newsome et al. approved that Fe nanoparticles were highly active in the recycling of technetium from groundwater in sediment systems [5]. Yazdani et al. demonstrated that the synthesis of Fe nanoparticle could be achieved by the drinking water treatment sludge method [6]. Eremin et al. examined another production approach by laser photolysis of Fe(CO) 5 vapor in the environment of argon, besides, the evaporation temperature of Fe nanoparticle is measured to be around 2100-2600 K [7]. From the commercial point of view, only the development of production technique is far not enough. This is because Fe and its alloys' products are prone to be corroded in nature. It was reported that the corrosion of metal surface has led huge loss of economy, and the efficient resistance of corrosion effects at natural environment is still a challenging task [8]. Obviously, the demand of surface protection for Fe nanoparticles is even larger than other macroscopic products due to their high surface-to-volume ratio.
Previously, the coating technology of Fe nanoparticles has attracted many researchers. Leng et al. dispersed Fe nanoparticles by polyvinyl pyrrolidone and found that silica layers could be directly coated without assistance from other surfactants [9]. While Dulinska-Molak et al. are more targeted on the application of biomedicines. They reported a successful synthesize and characterization of carbon-encapsulated Fe nanoparticles [10]. Khramtsov et al. further developed carbon coated Fe nanoparticles to optimize its conjugation with streptavidin and monoclonal antibodies [11]. Apart from inorganic coatings, organic polymer coatings are feasible for Fe nanoparticles as well. In order to keep highly-effective magnetic properties of Fe nanoparticles, McGrath et al. performed coating treatments basing on custom-synthesized phosphonate-grafted polyelectrolytes with different chain length [12]. Besides, other available polymer materials also include polyethylene glycol and polytetrahydrofuran. Galdames et al. generally applied Fe nanoparticles with those two kinds of coatings, respectively, for the application of lindane degradation.
It seems that previous coating studies for Fe nanoparticles mainly depended on experimental methods. Meanwhile, numerical simulations also play an irreplaceable role in surface coating or adsorption investigations. For example, MD simulations have been applied to model coating shell materials including carbon [13,14], nickel [15], ether [16], and ethanol [17,18]. Additionally, density functional theory (DFT) is a powerful tool to study the metal surface as well. Adsorbates like hydrogen [19,20], ethanol [21], oxygen [22,23], phosphorus, sulphur, and nitrogen [24] have been studied by DFT previously. The limitation of conventional DFT method is that issues like long-rang interactions were missed in its results. In this case, some studies aimed to make corrections by considering long-range effects [25][26][27].
Here we performed both MD simulations and DFT calculations to investigate the formation of ether coating layer on Fe nanoparticle surface. It was expected that those ether coated Fe nanoparticles would remain high activity and keep their qualified performances in conventional applications of Fe nanoparticles. In particular, MD simulations aimed to simulate the overall coating process from beginning to the end, and then formed a stably coated novel configuration. Next, DFT calculations targeted the optimized adsorption site and energy in this situation. Results of this study will provide basial support for industrial generations. It is believed that the formation of ether coating layer will effectively isolate Fe surface from atmosphere, and not affect industrial performances of Fe nanoparticles. The passivation effect of ether layer will effectively enlarge the retention period of Fe nanoparticles, because the corrosion effect for Fe nanoparticles is hindered.

MD Simulation and ReaxFF Force Field
With excellent calculation capabilities of microscope scale, MD simulations can effectively describe intermolecular behaviors by tracking their time dependent trajectories. The solution of MD simulations depends on numerical integrations of Newton's classical motion equations. So that movements, like acceleration, are determined by forces acting on those atoms or molecules. Specifically, MD simulations showed their highlights on obtaining structures, properties, and dynamic behaviors of nano-systems [28,29]. Updates of dynamic information including position, velocity and accelerations are proceeded and varied with time during the overall simulations [30]. In the present study, the calculation and integration of equilibrated MD simulations are performed by using Large Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, Sandia National Laboratories, Albuquerque, NX, USA) code package [31]. In which, the integration of motion equations is computed through Velocity-Verlet algorithm [32]. It was reported that LAMMPS showed its accuracy and efficiency in MD simulations of complex chemical systems [33].
Other than classical mechanics, force field in MD simulations refers to the potential functions used to describe interatomic interactions, such as potential energy and chemical bonds. The concept of bond order was introduced by Tersoff [34], and the determination of bond order is computed from the update of interatomic spaces at each timestep [18]. This allows continuous changes of angle, torsion, and atomic charge alongside with the break and form of chemical bonds [35]. Therefore, ReaxFF based MD simulation is a powerful instrument for heterogeneous systems [18]. In addition, the ReaxFF has also considered effects of polarization, such as the electronegativity equalization method, has been applied for geometry-dependent charge distributions [36]. Generally, ReaxFF obtains a system's total energy by combining all energy terms, as shown in Equation (1) [37].
where E system is the total energy; E bond is the bond energy; E over is the overcoordination energy; E under is the undercoordination energy; E lp is the lone-pair energy; E val is the valence angle energy; E tors is the torsion angle energy; E vdWaals is the van der Waals energy; and E coulomb is the coulomb energy. In this study, the specified ReaxFF force field was provided by Aryanpour et al. It was reported that the force field was initially developed for Fe-oxyhydroxide systems, and the training of its parameters has been validated by reproducing chemical reactions which showed good agreement with quantum mechanism and experimental investigations. Supplementary details of this ReaxFF relating to its training and setting can be obtained from Aryanpour, van Duin, and their co-workers [38].

MD Simulation Systems and Setup
Due to the demand of modelling and simulating surface coating behaviors for Fe nanoparticles, a standard Fe nanosphere was built initially. In terms of the equilibrium Fe lattice, it applied the body centered cubic (bcc) with lattice parameter set as 0.2866 nm [23]. Statistically, this sphere consisted of 2863 Fe atoms and its volume was about 12.56 nm 3 . Then the obtained model was placed at the center of a cube measuring 9 nm × 9 nm × 9 nm. As guided by relative simulations before [39][40][41], periodic boundaries were set along 3-Dimensional edges of the cube. Nevertheless, this initial model could be considered as a sphere whittled from Fe bulk crystal, and which might bring some unexpected edge effects. Therefore, the annealing treatment turned out to be necessary for this study. Generally, the change of temperature during the annealing simulation is presented in Figure 1 below. In the first stage of annealing process, the Fe nanoparticle was heated rapidly at a rate of 10.75 K/ps. The peak of heating temperature was 2150 K, and this temperature was already higher than experimental and computational melting points of bulk Fe [42]. Next, a freezing process was induced until the temperature of Fe nanoparticle was as low as 1700 K. At last, the freezing process was substituted by a slow cooling process. The cooling rate was about 1.65 K/ps, and the process ended when its temperature dropped to 0.1 K. Other than classical mechanics, force field in MD simulations refers to the potential functions used to describe interatomic interactions, such as potential energy and chemical bonds. The concept of bond order was introduced by Tersoff [34], and the determination of bond order is computed from the update of interatomic spaces at each timestep [18]. This allows continuous changes of angle, torsion, and atomic charge alongside with the break and form of chemical bonds [35]. Therefore, ReaxFF based MD simulation is a powerful instrument for heterogeneous systems [18]. In addition, the ReaxFF has also considered effects of polarization, such as the electronegativity equalization method, has been applied for geometry-dependent charge distributions [36]. Generally, ReaxFF obtains a system's total energy by combining all energy terms, as shown in Equation (1) [37]. (1) where Esystem is the total energy; Ebond is the bond energy; Eover is the overcoordination energy; Eunder is the undercoordination energy; Elp is the lone-pair energy; Eval is the valence angle energy; Etors is the torsion angle energy; EvdWaals is the van der Waals energy; and Ecoulomb is the coulomb energy. In this study, the specified ReaxFF force field was provided by Aryanpour et al. It was reported that the force field was initially developed for Fe-oxyhydroxide systems, and the training of its parameters has been validated by reproducing chemical reactions which showed good agreement with quantum mechanism and experimental investigations. Supplementary details of this ReaxFF relating to its training and setting can be obtained from Aryanpour, van Duin, and their co-workers [38].

MD Simulation Systems and Setup
Due to the demand of modelling and simulating surface coating behaviors for Fe nanoparticles, a standard Fe nanosphere was built initially. In terms of the equilibrium Fe lattice, it applied the body centered cubic (bcc) with lattice parameter set as 0.2866 nm [23]. Statistically, this sphere consisted of 2863 Fe atoms and its volume was about 12.56 nm 3 . Then the obtained model was placed at the center of a cube measuring 9 nm × 9 nm × 9 nm. As guided by relative simulations before [39][40][41], periodic boundaries were set along 3-Dimensional edges of the cube. Nevertheless, this initial model could be considered as a sphere whittled from Fe bulk crystal, and which might bring some unexpected edge effects. Therefore, the annealing treatment turned out to be necessary for this study. Generally, the change of temperature during the annealing simulation is presented in Figure 1 below. In the first stage of annealing process, the Fe nanoparticle was heated rapidly at a rate of 10.75 K/ps. The peak of heating temperature was 2150 K, and this temperature was already higher than experimental and computational melting points of bulk Fe [42]. Next, a freezing process was induced until the temperature of Fe nanoparticle was as low as 1700 K. At last, the freezing process was substituted by a slow cooling process. The cooling rate was about 1.65 K/ps, and the process ended when its temperature dropped to 0.1 K.  As one of the expected innovations in this study, a novel double-sphere model was induced to prove whether coating behaviors was affected by the gap between two Fe nanoparticles or not. In this study, two annealed Fe nanoparticles were paced together in a 9 nm × 13 nm × 9 nm scaled simulation box. The horizontal gap distance between those two spheres was equal to 0 nm, in other words, they were just connected with each other. As for the settlement of ether molecules, all of them have been geometrically optimized by DFT methods (discussed in Section 3.1). Then, the packing optimization of initial liquid ether configuration was achieved by the Packmol software (version 18.169) package [43]. According to the real bulk density of liquid ether at room temperature, there were 2908 ether molecules set for the single sphere model and 3287 ether molecules for double-sphere model. Next, Fe nanoparticles and ether molecules were assembled by visual molecular dynamics (VMD) [44]. Due to the consideration of coating speed and effect [17,39], the temperature of those ether molecules were kept at about 750 K during the whole simulation. In MD simulations, the temperature of model can be modified by at least three theories, such as Nose [45]/Hoover [46], Berendsen [47], and Anderson [48]. It was reported that the canonical ensemble (NVT) with Nose/Hoover thermostat method can truly obtain canonical distribution of particles' dynamics [40]. Therefore, all temperatures in this study were controlled by NVT theory with a damping parameter of 25 fs. For both of those ether coating systems, a total of 4 × 10 6 integrations, each at 0.25 fs, were executed [38]. In addition, all MD configurations of studied models were rendered by the OVITO software package (version: 2.9.0; http://ovito.org/) [49].

First-Principle Calculations
Benefited by the development of computational capabilities, DFT method can provide accurate and robust assistance for nanoscale research and applications [50]. In this study, spin-dependent DFT calculations have been implemented by PWscf code of the Quantum Espresso (QE) program [51,52], and DFT configurations are presented via VESTA [53]. Ultrasoft pseudopotentials were provided by the official website of QE, and they were performed on the Perdew-Burke-Ernzerhof (GGA-PBE) level of approximation [54]. The cutoff kinetic energy of plane waves was valued as 30 Ry, and which was about 1/10 of charge density cutoff [19]. Through Monckhorst-Pack scheme, a (7 × 7 × 1) series of K points was set for Brillouin zone. In contrast with MD simulations, the equilibrium lattice constant used in DFT has been recalculated as about 2.83 Å, which was in good agreement with other reported results [23]. Then the value was applied to building the (111) type of Fe surface with seven layers. Along the vertical direction of the cell, there was a vacuum space with 16 Å high separating periodical Fe (110) slabs. As shown in Figure 2, the adsorption of Fe (110) surface was studied at four different adsorption sites, including threefold site (hcp), twofold site (long bridge and short bridge site), and onefold site (atop). The coverage ratio in this case is about 0.6 ML.
Coatings 2019, 9, x FOR PEER REVIEW 4 of 14 As one of the expected innovations in this study, a novel double-sphere model was induced to prove whether coating behaviors was affected by the gap between two Fe nanoparticles or not. In this study, two annealed Fe nanoparticles were paced together in a 9 nm × 13 nm × 9 nm scaled simulation box. The horizontal gap distance between those two spheres was equal to 0 nm, in other words, they were just connected with each other. As for the settlement of ether molecules, all of them have been geometrically optimized by DFT methods (discussed in Section 3.1). Then, the packing optimization of initial liquid ether configuration was achieved by the Packmol software (version 18.169) package [43]. According to the real bulk density of liquid ether at room temperature, there were 2908 ether molecules set for the single sphere model and 3287 ether molecules for double-sphere model. Next, Fe nanoparticles and ether molecules were assembled by visual molecular dynamics (VMD) [44]. Due to the consideration of coating speed and effect [17,39], the temperature of those ether molecules were kept at about 750 K during the whole simulation. In MD simulations, the temperature of model can be modified by at least three theories, such as Nose [45]/Hoover [46], Berendsen [47], and Anderson [48]. It was reported that the canonical ensemble (NVT) with Nose/Hoover thermostat method can truly obtain canonical distribution of particles' dynamics [40]. Therefore, all temperatures in this study were controlled by NVT theory with a damping parameter of 25 fs. For both of those ether coating systems, a total of 4 × 10 6 integrations, each at 0.25 fs, were executed [38]. In addition, all MD configurations of studied models were rendered by the OVITO software package (version: 2.9.0; http://ovito.org/) [49].

First-Principle Calculations.
Benefited by the development of computational capabilities, DFT method can provide accurate and robust assistance for nanoscale research and applications [50]. In this study, spin-dependent DFT calculations have been implemented by PWscf code of the Quantum Espresso (QE) program [51,52], and DFT configurations are presented via VESTA [53]. Ultrasoft pseudopotentials were provided by the official website of QE, and they were performed on the Perdew-Burke-Ernzerhof (GGA-PBE) level of approximation [54]. The cutoff kinetic energy of plane waves was valued as 30 Ry, and which was about 1/10 of charge density cutoff [19]. Through Monckhorst-Pack scheme, a (7 × 7 × 1) series of K points was set for Brillouin zone. In contrast with MD simulations, the equilibrium lattice constant used in DFT has been recalculated as about 2.83 Å, which was in good agreement with other reported results [23]. Then the value was applied to building the (111) type of Fe surface with seven layers. Along the vertical direction of the cell, there was a vacuum space with 16 Å high separating periodical Fe (110) slabs. As shown in Figure 2, the adsorption of Fe (110) surface was studied at four different adsorption sites, including threefold site (hcp), twofold site (long bridge and short bridge site), and onefold site (atop). The coverage ratio in this case is about 0.6 ML. For both metal surface and ether molecule, their relaxations were considered to be converged while forces acting on them are not higher than 0.05 eV/Å. Due to the demand of stability analysis for adsorbed configurations, Equation (2) is used to compute the adsorption energy of each model.
(2) For both metal surface and ether molecule, their relaxations were considered to be converged while forces acting on them are not higher than 0.05 eV/Å. Due to the demand of stability analysis for adsorbed configurations, Equation (2) is used to compute the adsorption energy of each model.
where E ads is the expected adsorption energy of each adsorbing model; E Fe + Ether , E Fe , and E Ether are total energies of assembled model, Fe (110) supercell, and single ether molecule, respectively.

First-Principle Calculations
As the adsorbate in this study, geometrical optimization and energetical minimization of ether molecule play an elemental role for the following DFT and MD simulations. Figure 3 presents the chemical structure of ether molecule, in which all atoms of ether molecule have been labelled, respectively, to associate each atom with its charge value. According to the DFT population analysis result, the charge of each atom is listed in Table 1. Later, those parameters will be set as initial datafile for MD simulations. Next, an optimized ether molecule is placed above the Fe (110) surface to obtain the optimal adsorption site. As mentioned above, there are four potential sites applied to calculate the adsorption energy, and the initial distance between O1 atom and the first layer of Fe (110) is about 3.0 Å. Taking the long bridge site as an example, the adsorbed configuration is presented as Figure 4. Through DFT calculation, the adsorption energy for each site is summarized in Table 2. As shown in the table, all those four computed adsorption energies are negative values meaning that this is an exothermic adsorption process. Thus, it could be deduced that more energy will be released if the adsorption is stronger. According to their comparisons, it appears that the long bridge adsorption site on Fe (110) is an optimal site for the adsorption of ether molecule. In other words, the obtained configuration will be more stable than other configurations if an ether molecule is adsorbed on long bridge site. Finally, details of original and adsorbed ether molecules' chemical bonds are summarised in Table 3. where Eads is the expected adsorption energy of each adsorbing model; EFe + Ether, EFe, and EEther are total energies of assembled model, Fe (110) supercell, and single ether molecule, respectively.

First-Principle Calculations
As the adsorbate in this study, geometrical optimization and energetical minimization of ether molecule play an elemental role for the following DFT and MD simulations. Figure 3 presents the chemical structure of ether molecule, in which all atoms of ether molecule have been labelled, respectively, to associate each atom with its charge value. According to the DFT population analysis result, the charge of each atom is listed in Table 1. Later, those parameters will be set as initial datafile for MD simulations. Next, an optimized ether molecule is placed above the Fe (110) surface to obtain the optimal adsorption site. As mentioned above, there are four potential sites applied to calculate the adsorption energy, and the initial distance between O1 atom and the first layer of Fe (110) is about 3.0 Å. Taking the long bridge site as an example, the adsorbed configuration is presented as Figure 4. Through DFT calculation, the adsorption energy for each site is summarized in Table 2. As shown in the table, all those four computed adsorption energies are negative values meaning that this is an exothermic adsorption process. Thus, it could be deduced that more energy will be released if the adsorption is stronger. According to their comparisons, it appears that the long bridge adsorption site on Fe (110) is an optimal site for the adsorption of ether molecule. In other words, the obtained configuration will be more stable than other configurations if an ether molecule is adsorbed on long bridge site. Finally, details of original and adsorbed ether molecules' chemical bonds are summarised in Table 3.     where Eads is the expected adsorption energy of each adsorbing model; EFe + Ether, EFe, and EEther are total energies of assembled model, Fe (110) supercell, and single ether molecule, respectively.

First-Principle Calculations
As the adsorbate in this study, geometrical optimization and energetical minimization of ether molecule play an elemental role for the following DFT and MD simulations. Figure 3 presents the chemical structure of ether molecule, in which all atoms of ether molecule have been labelled, respectively, to associate each atom with its charge value. According to the DFT population analysis result, the charge of each atom is listed in Table 1. Later, those parameters will be set as initial datafile for MD simulations. Next, an optimized ether molecule is placed above the Fe (110) surface to obtain the optimal adsorption site. As mentioned above, there are four potential sites applied to calculate the adsorption energy, and the initial distance between O1 atom and the first layer of Fe (110) is about 3.0 Å. Taking the long bridge site as an example, the adsorbed configuration is presented as Figure 4. Through DFT calculation, the adsorption energy for each site is summarized in Table 2. As shown in the table, all those four computed adsorption energies are negative values meaning that this is an exothermic adsorption process. Thus, it could be deduced that more energy will be released if the adsorption is stronger. According to their comparisons, it appears that the long bridge adsorption site on Fe (110) is an optimal site for the adsorption of ether molecule. In other words, the obtained configuration will be more stable than other configurations if an ether molecule is adsorbed on long bridge site. Finally, details of original and adsorbed ether molecules' chemical bonds are summarised in Table 3.

Single Sphere Coating Model
To uncover the mechanism of coating Fe nanoparticles by liquid ether, a single sphere model has been placed to be surrounded by the ether molecules. Although the model was performed for totally 1 ns, it is observed that the coating behavior mainly occurred in the first 30 ps. Figure 5 exhibits four snapshots of adsorption configurations with intervals of 8 ps. Those snapshots are cross-sections of the single sphere model and each is about 25 Å thick. Initially (0 ps), the Fe nanoparticle and bulk of ether molecules were separated by a vacuum about 15 Å long. After 8 ps, several ether molecules were moved into the vacuum and captured on Fe surface. During this period of time, the density of bulk region was still higher than vacuum. In the next two intervals, the distribution of ether molecules tended to be uniform. While, in this process, most parts of Fe surface were coated successfully. However, the coating effect was not isotropy, especially because the bottom right corner of Fe particle was explicitly coated less than other parts. Therefore, it could be summarized that the coating behavior on this occasion did not occur simultaneously in 360 • , but spread step by step until the distribution was stabilized. The stability of adsorption system was not convincing enough to be judged visually. Therefore, a statistic of potential energy within 200 ps is plotted, as shown in Figure 6. The figure illustrates that the potential energy of system decreased drastically during the early stage of adsorption, then stabilized from 45 ps. Also after 150 ps, the system was reasonable to be considered as an equilibrated system.

Single Sphere Coating Model
To uncover the mechanism of coating Fe nanoparticles by liquid ether, a single sphere model has been placed to be surrounded by the ether molecules. Although the model was performed for totally 1 ns, it is observed that the coating behavior mainly occurred in the first 30 ps. Figure 5 exhibits four snapshots of adsorption configurations with intervals of 8 ps. Those snapshots are cross-sections of the single sphere model and each is about 25 Å thick. Initially (0 ps), the Fe nanoparticle and bulk of ether molecules were separated by a vacuum about 15 Å long. After 8 ps, several ether molecules were moved into the vacuum and captured on Fe surface. During this period of time, the density of bulk region was still higher than vacuum. In the next two intervals, the distribution of ether molecules tended to be uniform. While, in this process, most parts of Fe surface were coated successfully. However, the coating effect was not isotropy, especially because the bottom right corner of Fe particle was explicitly coated less than other parts. Therefore, it could be summarized that the coating behavior on this occasion did not occur simultaneously in 360°, but spread step by step until the distribution was stabilized. The stability of adsorption system was not convincing enough to be judged visually. Therefore, a statistic of potential energy within 200 ps is plotted, as shown in Figure  6. The figure illustrates that the potential energy of system decreased drastically during the early stage of adsorption, then stabilized from 45 ps. Also after 150 ps, the system was reasonable to be considered as an equilibrated system.   Despite the plot of potential energy, the adsorption curve can also be used as a tool to analyze the adsorption process in detail. As Figure 7 displays, the adsorption of ether molecules could be basically divided into three stages. Stage 1 refers to the first 30 ps, in which ether molecules were adsorbed effectively and the curve sharply climbed up to the quantity of 280 molecules. Then, from 30 to 45 ps, the rate of adsorption slowed down, and the curve became smoother. Finally, the stabilized adsorption curve remains fluctuating around the quantity of 340 molecules. Figure 8 depicts the charge distribution of Fe nanoparticle at 1 ns when all ether molecules have been removed. Similarly, ¼ of particle have been sliced with the aim to present the internal distribution. It confirms that the adsorption process was typically a charge transfer process as well. Charge values of Fe atoms were almost zero initially, and then negative charges moved from Fe to ether atoms with the occurrence of adsorption. Therefore, the transfer and attraction effect between Fe atoms and ether molecules generally motivated adsorption behaviors.   Despite the plot of potential energy, the adsorption curve can also be used as a tool to analyze the adsorption process in detail. As Figure 7 displays, the adsorption of ether molecules could be basically divided into three stages. Stage 1 refers to the first 30 ps, in which ether molecules were adsorbed effectively and the curve sharply climbed up to the quantity of 280 molecules. Then, from 30 to 45 ps, the rate of adsorption slowed down, and the curve became smoother. Finally, the stabilized adsorption curve remains fluctuating around the quantity of 340 molecules. Figure 8 depicts the charge distribution of Fe nanoparticle at 1 ns when all ether molecules have been removed. Similarly, 1 4 of particle have been sliced with the aim to present the internal distribution. It confirms that the adsorption process was typically a charge transfer process as well. Charge values of Fe atoms were almost zero initially, and then negative charges moved from Fe to ether atoms with the occurrence of adsorption. Therefore, the transfer and attraction effect between Fe atoms and ether molecules generally motivated adsorption behaviors. Despite the plot of potential energy, the adsorption curve can also be used as a tool to analyze the adsorption process in detail. As Figure 7 displays, the adsorption of ether molecules could be basically divided into three stages. Stage 1 refers to the first 30 ps, in which ether molecules were adsorbed effectively and the curve sharply climbed up to the quantity of 280 molecules. Then, from 30 to 45 ps, the rate of adsorption slowed down, and the curve became smoother. Finally, the stabilized adsorption curve remains fluctuating around the quantity of 340 molecules. Figure 8 depicts the charge distribution of Fe nanoparticle at 1 ns when all ether molecules have been removed. Similarly, ¼ of particle have been sliced with the aim to present the internal distribution. It confirms that the adsorption process was typically a charge transfer process as well. Charge values of Fe atoms were almost zero initially, and then negative charges moved from Fe to ether atoms with the occurrence of adsorption. Therefore, the transfer and attraction effect between Fe atoms and ether molecules generally motivated adsorption behaviors.    Despite the plot of potential energy, the adsorption curve can also be used as a tool to analyze the adsorption process in detail. As Figure 7 displays, the adsorption of ether molecules could be basically divided into three stages. Stage 1 refers to the first 30 ps, in which ether molecules were adsorbed effectively and the curve sharply climbed up to the quantity of 280 molecules. Then, from 30 to 45 ps, the rate of adsorption slowed down, and the curve became smoother. Finally, the stabilized adsorption curve remains fluctuating around the quantity of 340 molecules. Figure 8 depicts the charge distribution of Fe nanoparticle at 1 ns when all ether molecules have been removed. Similarly, ¼ of particle have been sliced with the aim to present the internal distribution. It confirms that the adsorption process was typically a charge transfer process as well. Charge values of Fe atoms were almost zero initially, and then negative charges moved from Fe to ether atoms with the occurrence of adsorption. Therefore, the transfer and attraction effect between Fe atoms and ether molecules generally motivated adsorption behaviors.   Besides, such a strong attraction effect from plenty of ether molecules with high kinetic energy will inevitably cause radial pressure on Fe nanoparticle. Figure 9 is the analysis of Fe displacement at 200 ps. According to this figure, it is observed that pressure from the bulk of ether made Fe atoms left their original coordinates. Typically, color coded configuration shows that shell atoms were explicitly affected larger than core atoms. Also, there was little displacement within the radius of 18 Å. The closer the Fe atoms were to shell surface, the higher the extent of displacement would be. When the maximum displacement reached for about 6.36 Å, the radius of particle would be enlarged until 25.66 Å. Besides, Figure 10 displaces the dynamic distribution of shell atoms' displacement. From 0 to 200 ps, those displacement behaviors were mainly happened before 45 ps. The most important aspect is that this result is in line with the adsorption cure ( Figure 7). As such, the adsorption of ether molecules is confirmed to effectively cause displacements of shell Fe atoms. Besides, such a strong attraction effect from plenty of ether molecules with high kinetic energy will inevitably cause radial pressure on Fe nanoparticle. Figure 9 is the analysis of Fe displacement at 200 ps. According to this figure, it is observed that pressure from the bulk of ether made Fe atoms left their original coordinates. Typically, color coded configuration shows that shell atoms were explicitly affected larger than core atoms. Also, there was little displacement within the radius of 18 Å. The closer the Fe atoms were to shell surface, the higher the extent of displacement would be. When the maximum displacement reached for about 6.36 Å, the radius of particle would be enlarged until 25.66 Å. Besides, Figure 10 displaces the dynamic distribution of shell atoms' displacement. From 0 to 200 ps, those displacement behaviors were mainly happened before 45 ps. The most important aspect is that this result is in line with the adsorption cure ( Figure 7). As such, the adsorption of ether molecules is confirmed to effectively cause displacements of shell Fe atoms.  Until then, the preparation of ether coated Fe nanoparticle is still unfinished, especially since the Fe nanoparticle is still soaked under liquid ether. The final step will be a filtering treatment in order to form a purely coated particle without other free molecules. Previously, we successfully filtered ether coated Al nanoparticles [16]. This time, the Fe nanoparticle system was treated with the same principle and method. Figure 11 presents the desorption process caused by filtering treatment. The figure demonstrates that the ether coated configuration has been stabilized at room temperature (300 K), because there was no desorbed ether molecule since cycle 8. Figure 11 also displays the obtained coating configuration with 259 ether molecules. Besides, such a strong attraction effect from plenty of ether molecules with high kinetic energy will inevitably cause radial pressure on Fe nanoparticle. Figure 9 is the analysis of Fe displacement at 200 ps. According to this figure, it is observed that pressure from the bulk of ether made Fe atoms left their original coordinates. Typically, color coded configuration shows that shell atoms were explicitly affected larger than core atoms. Also, there was little displacement within the radius of 18 Å. The closer the Fe atoms were to shell surface, the higher the extent of displacement would be. When the maximum displacement reached for about 6.36 Å, the radius of particle would be enlarged until 25.66 Å. Besides, Figure 10 displaces the dynamic distribution of shell atoms' displacement. From 0 to 200 ps, those displacement behaviors were mainly happened before 45 ps. The most important aspect is that this result is in line with the adsorption cure ( Figure 7). As such, the adsorption of ether molecules is confirmed to effectively cause displacements of shell Fe atoms.  Until then, the preparation of ether coated Fe nanoparticle is still unfinished, especially since the Fe nanoparticle is still soaked under liquid ether. The final step will be a filtering treatment in order to form a purely coated particle without other free molecules. Previously, we successfully filtered ether coated Al nanoparticles [16]. This time, the Fe nanoparticle system was treated with the same principle and method. Figure 11 presents the desorption process caused by filtering treatment. The figure demonstrates that the ether coated configuration has been stabilized at room temperature (300 K), because there was no desorbed ether molecule since cycle 8. Figure 11 also displays the obtained coating configuration with 259 ether molecules. Until then, the preparation of ether coated Fe nanoparticle is still unfinished, especially since the Fe nanoparticle is still soaked under liquid ether. The final step will be a filtering treatment in order to form a purely coated particle without other free molecules. Previously, we successfully filtered ether coated Al nanoparticles [16]. This time, the Fe nanoparticle system was treated with the same principle and method. Figure 11 presents the desorption process caused by filtering treatment. The figure demonstrates that the ether coated configuration has been stabilized at room temperature (300 K), because there was no desorbed ether molecule since cycle 8. Figure 11 also displays the obtained coating configuration with 259 ether molecules.

Double Sphere Coating Model
Although the single sphere model has been broadly applied to study dynamic or reaction behaviors of various metal nanoparticles previously [14,17,55], practical impacts from multiparticle environment is still inevitable. In this study, a double sphere model was built to further uncover coating mechanisms for multiparticle model. As mentioned previously, double Fe particles model have been pre-built, and the settlement of ether molecules was entirely same with single sphere model (Section 3.2). Figure 12 is the cross-section snapshot of double sphere model from 0 to 15 ps. Similar with the single sphere model, the coating behavior in this section was neither performed by uniform distribution. At about 4 ps, there were two vacuum spaces formed symmetrically between two nanoparticles. Then their volume shrank step by step and were finally filled up after 8 ps. According to the configuration at 15 ps, all parts of Fe shell have been coated entirely and the appearance of neighbor particles did not lead to exposure of particle surface. As the evaluation of coating process, another adsorption curve is plotted as shown in Figure 13. In which, stage 1 refers to the rapid coating process happened from 0 to 15 ps, then the coating curve slowly trended to be stabilized until 40 ps. Generally, the time length of overall coating process was

Double Sphere Coating Model
Although the single sphere model has been broadly applied to study dynamic or reaction behaviors of various metal nanoparticles previously [14,17,55], practical impacts from multiparticle environment is still inevitable. In this study, a double sphere model was built to further uncover coating mechanisms for multiparticle model. As mentioned previously, double Fe particles model have been pre-built, and the settlement of ether molecules was entirely same with single sphere model (Section 3.2). Figure 12 is the cross-section snapshot of double sphere model from 0 to 15 ps. Similar with the single sphere model, the coating behavior in this section was neither performed by uniform distribution. At about 4 ps, there were two vacuum spaces formed symmetrically between two nanoparticles. Then their volume shrank step by step and were finally filled up after 8 ps. According to the configuration at 15 ps, all parts of Fe shell have been coated entirely and the appearance of neighbor particles did not lead to exposure of particle surface.

Double Sphere Coating Model
Although the single sphere model has been broadly applied to study dynamic or reaction behaviors of various metal nanoparticles previously [14,17,55], practical impacts from multiparticle environment is still inevitable. In this study, a double sphere model was built to further uncover coating mechanisms for multiparticle model. As mentioned previously, double Fe particles model have been pre-built, and the settlement of ether molecules was entirely same with single sphere model (Section 3.2). Figure 12 is the cross-section snapshot of double sphere model from 0 to 15 ps. Similar with the single sphere model, the coating behavior in this section was neither performed by uniform distribution. At about 4 ps, there were two vacuum spaces formed symmetrically between two nanoparticles. Then their volume shrank step by step and were finally filled up after 8 ps. According to the configuration at 15 ps, all parts of Fe shell have been coated entirely and the appearance of neighbor particles did not lead to exposure of particle surface. As the evaluation of coating process, another adsorption curve is plotted as shown in Figure 13. In which, stage 1 refers to the rapid coating process happened from 0 to 15 ps, then the coating curve slowly trended to be stabilized until 40 ps. Generally, the time length of overall coating process was As the evaluation of coating process, another adsorption curve is plotted as shown in Figure 13. In which, stage 1 refers to the rapid coating process happened from 0 to 15 ps, then the coating curve slowly trended to be stabilized until 40 ps. Generally, the time length of overall coating process was in line with the single sphere model. From 40 ps to 1 ns, there were totally about 638 ether molecules coated on Fe surface. The proportion of maximum adsorptions for double sphere model and single sphere model was about 1.876. Due to the consideration of overlapped area, this proportion means that the coating effect for those two models are similar with each other. Analyzed from the angle of potential energy, the equilibrium state of double sphere system is represented by Figure 14. To further investigate the impact of multiparticle environment, uncontacted and contacted regions of the left nanoparticle are labelled as regions A and B respectively as indicated in Figure 15. Then, each of their adsorption curves is plotted as well. It is observed that the adsorption behaviour of ether molecules for region B was typically affected. This is because parts of particle surface are covered by the other one, and then the availability of adsorption site is limited. Due to the consideration of overlapped area, this proportion means that the coating effect for those two models are similar with each other. Analyzed from the angle of potential energy, the equilibrium state of double sphere system is represented by Figure 14. To further investigate the impact of multiparticle environment, uncontacted and contacted regions of the left nanoparticle are labelled as regions A and B respectively as indicated in Figure 15. Then, each of their adsorption curves is plotted as well. It is observed that the adsorption behaviour of ether molecules for region B was typically affected. This is because parts of particle surface are covered by the other one, and then the availability of adsorption site is limited.    Due to the consideration of overlapped area, this proportion means that the coating effect for those two models are similar with each other. Analyzed from the angle of potential energy, the equilibrium state of double sphere system is represented by Figure 14. To further investigate the impact of multiparticle environment, uncontacted and contacted regions of the left nanoparticle are labelled as regions A and B respectively as indicated in Figure 15. Then, each of their adsorption curves is plotted as well. It is observed that the adsorption behaviour of ether molecules for region B was typically affected. This is because parts of particle surface are covered by the other one, and then the availability of adsorption site is limited.     The proportion of maximum adsorptions for double sphere model and single sphere model was about 1.876. Due to the consideration of overlapped area, this proportion means that the coating effect for those two models are similar with each other. Analyzed from the angle of potential energy, the equilibrium state of double sphere system is represented by Figure 14. To further investigate the impact of multiparticle environment, uncontacted and contacted regions of the left nanoparticle are labelled as regions A and B respectively as indicated in Figure 15. Then, each of their adsorption curves is plotted as well. It is observed that the adsorption behaviour of ether molecules for region B was typically affected. This is because parts of particle surface are covered by the other one, and then the availability of adsorption site is limited.     Figure 16 is the cross section of double sphere model at 1 ns. The figure is colored by the extent of the atoms' movement in order to present their displacement distribution. It could be obviously observed that a layer of ether molecules with relatively low displacement was coated around the shell of Fe nanoparticles. In other words, those molecules were proved to be captured by the adsorption effects. In contrast with the single sphere model, the gap between two particles was under the double adsorption influence. Dynamic properties of ether molecules in this region ought to be different from other free molecules, although they were not entirely adsorbed. Simultaneously affected by both of those nanoparticles, the region between two nanoparticles were more likely to be with low displacement and the trajectory of their movement was restricted in this situation. Kinetic energy is another property to evaluate the double adsorption influence. This is because the value of kinetic energy is directly determined by the velocity of each atom. Just as in Figure 16, the distribution of kinetic energy at the same moment is presented in Figure 17. In this figure, the double adsorption influence is presented more clearly than in Figure 16, because the low kinetic energy area between two nanoparticles is explicitly thicker than other parts of configuration.
Coatings 2019, 9, x FOR PEER REVIEW 11 of 14 Figure 16 is the cross section of double sphere model at 1 ns. The figure is colored by the extent of the atoms' movement in order to present their displacement distribution. It could be obviously observed that a layer of ether molecules with relatively low displacement was coated around the shell of Fe nanoparticles. In other words, those molecules were proved to be captured by the adsorption effects. In contrast with the single sphere model, the gap between two particles was under the double adsorption influence. Dynamic properties of ether molecules in this region ought to be different from other free molecules, although they were not entirely adsorbed. Simultaneously affected by both of those nanoparticles, the region between two nanoparticles were more likely to be with low displacement and the trajectory of their movement was restricted in this situation. Kinetic energy is another property to evaluate the double adsorption influence. This is because the value of kinetic energy is directly determined by the velocity of each atom. Just as in Figure 16, the distribution of kinetic energy at the same moment is presented in Figure 17. In this figure, the double adsorption influence is presented more clearly than in Figure 16, because the low kinetic energy area between two nanoparticles is explicitly thicker than other parts of configuration.

Conclusion
Based on first-principle calculations and MD simulations, the mechanism of coating ether on Fe nanoparticles has been discussed in this study. Above all, the obtained DFT results show that this coating behavior was generally an exothermic process, because the energy of adsorbed system was reduced comparing with the separated components. Among four computed adsorption sites, the long bridge site was confirmed to be the most stable site for ether adsorption. According to MD simulations of the single sphere model, the general coating process was found to be non-uniform, as some parts of the adsorption surface were coated more rapidly than other parts. However, all surfaces were stably coated in the end. The equilibrium state is validated by plots of both the potential energy and adsorption curve. It is also demonstrated that the outer shell of Fe atoms was affected significantly by the coating effects, exhibiting changes such as atomic charge transfer and   Figure 16 is the cross section of double sphere model at 1 ns. The figure is colored by the extent of the atoms' movement in order to present their displacement distribution. It could be obviously observed that a layer of ether molecules with relatively low displacement was coated around the shell of Fe nanoparticles. In other words, those molecules were proved to be captured by the adsorption effects. In contrast with the single sphere model, the gap between two particles was under the double adsorption influence. Dynamic properties of ether molecules in this region ought to be different from other free molecules, although they were not entirely adsorbed. Simultaneously affected by both of those nanoparticles, the region between two nanoparticles were more likely to be with low displacement and the trajectory of their movement was restricted in this situation. Kinetic energy is another property to evaluate the double adsorption influence. This is because the value of kinetic energy is directly determined by the velocity of each atom. Just as in Figure 16, the distribution of kinetic energy at the same moment is presented in Figure 17. In this figure, the double adsorption influence is presented more clearly than in Figure 16, because the low kinetic energy area between two nanoparticles is explicitly thicker than other parts of configuration.

Conclusion
Based on first-principle calculations and MD simulations, the mechanism of coating ether on Fe nanoparticles has been discussed in this study. Above all, the obtained DFT results show that this coating behavior was generally an exothermic process, because the energy of adsorbed system was reduced comparing with the separated components. Among four computed adsorption sites, the long bridge site was confirmed to be the most stable site for ether adsorption. According to MD simulations of the single sphere model, the general coating process was found to be non-uniform, as some parts of the adsorption surface were coated more rapidly than other parts. However, all surfaces were stably coated in the end. The equilibrium state is validated by plots of both the potential energy and adsorption curve. It is also demonstrated that the outer shell of Fe atoms was affected significantly by the coating effects, exhibiting changes such as atomic charge transfer and

Conclusions
Based on first-principle calculations and MD simulations, the mechanism of coating ether on Fe nanoparticles has been discussed in this study. Above all, the obtained DFT results show that this coating behavior was generally an exothermic process, because the energy of adsorbed system was reduced comparing with the separated components. Among four computed adsorption sites, the long bridge site was confirmed to be the most stable site for ether adsorption. According to MD simulations of the single sphere model, the general coating process was found to be non-uniform, as some parts of the adsorption surface were coated more rapidly than other parts. However, all surfaces were stably coated in the end. The equilibrium state is validated by plots of both the potential energy and adsorption curve. It is also demonstrated that the outer shell of Fe atoms was affected significantly by the coating effects, exhibiting changes such as atomic charge transfer and displacement. Finally, the purely coated configuration was successfully obtained. Furthermore, the double sphere model showed that the gap between two particles was coated later than other parts of particles. Besides, the number of coated ether molecules reduced as well. Meanwhile, the movement of uncoated ether molecules in this gap was hindered because of the double adsorption effects of two neighboring particles.