Computational Simulations of Nanoconﬁned Argon Film through Adsorption–Desorption in a Uniform Slit Pore

: We performed hybrid grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations to investigate the adsorption-desorption isotherms of argon molecules conﬁned between commensurate and incommensurate contacts in nanoscale thickness. The recently proposed mid-density scheme was applied to the obtained hysteresis loops to produce a realistic equilibrium phase of nanoconﬁned ﬂuids. The appropriate chemical potentials can be determined if the equilibrium structures predicted by GCMC/MD simulations are consistent with those observed in previously developed liquid-vapor molecular dynamics (LVMD) simulations. With the chemical potential as input, the equilibrium structures obtained by GCMC/MD simulations can be used as reasonable initial conﬁgurations for future metadynamics free energy calculations.


Introduction
The structure and thermodynamic state of a liquid film confined between two solid surfaces in nanoscale thickness and its mechanical response to external loading are very important topics in surface science [1], which has attracted considerable interests and research efforts for decades in many surface force experiments (mainly in surface force apparatus (SFA) or surface force balance (SFB) measurements) [2][3][4][5] and molecular simulations [6,7]. The fundamental questions concerning the squeezing and shearing behaviors of nanoconfined fluids are particularly relevant to the fields of nanotribology, biolubrication, coating, and surface engineering and have been a matter of persistent debate despite many new findings and new insights [8][9][10][11].
The very low compression rate or lateral sliding speed of mechanical driving springs in an SFA/B instrument largely enables a nanoconfined liquid film in a quasistatic equilibrium state, except in an unstable jump during the layering transition or in a fast slip during a stick-slip sliding. Because of this, we need to first understand the actual equilibrium structure of a nanoconfined liquid film that should correspond to a lowest free energy state. These equilibrium states may not be fully explored by the SFA/B mechanical driving system due to its nonequilibrium nature of driven dynamics during unstable transitions.
Recently, a liquid-vapor molecular dynamics (LVMD) ensemble was used to investigate the force oscillation, phase transition, and shearing behaviors of liquid films confined between two solid walls [12][13][14][15][16]. The main features of the LVMD simulation can be summarized in three aspects: First, it enables the squeeze-out of the confined fluid during normal compression due to the extended side walls along the squeeze-out direction. Second, the lateral pressure of the confined system is comparable to the ambient pressure because of the coexistence of vapor phase around the central liquid droplet. Third, the squeeze-out of the confined fluid is realized by compressing the driving spring, which is connected to the top confining wall. Such a mechanical loading is similar to the situation in SFA/B experiments. Therefore, an LVMD ensemble is particularly suitable in investigating the squeeze-out behavior of confined fluids in surface force experiments. It has also been extended to study the mechanical response of a liquid film confined between AFM tip and substrate [17][18][19]. A critical issue in the LVMD simulation of a nanoconfined fluid concerns whether the simulation can produce a realistic thermodynamic equilibrium state of the liquid film under a sufficiently low compression rate. In our previous study [14], grand canonical Monte Carlo (GCMC) and LVMD simulations were carried out to study the squeezing and phase transition of simple liquid argon confined between two commensurate solid surfaces. The comparative study revealed that LVMD simulation can capture a major part of the force oscillation profile predicted by the equilibrium GCMC simulation and, at the same time, reveal the non-equilibrium squeeze-out behavior of the confined film. Still, the problem of whether the confined liquid film is in the lowest free energy state remains.
Metadynamics [20,21] proves to be a state-of-the-art computational framework for enhanced sampling coupled with molecular dynamics (MD) simulations in an efficient, flexible, and accurate manner. Its ultimate goal is to simulate rare events and reconstruct the free energy surface (FES) as a function of collective variables (CVs) in physics/biophysics, chemistry, and material science [22,23]. To study phase transitions involving the solid nucleation and/or relative free energy of different phases, structural order parameters (such as Steinhardt parameters [24,25]) can be selected as CVs [26,27]. In this sense, however, LVMD assembly as the simulation setup is not suitable for a metadynamics study, due to the coexistence of multiple phases involved in LVMD simulation. The simulation setup in a uniform slit pore with periodic boundary conditions (PBCs) applied in the lateral directions (i.e., there is no liquid-gas or solid-liquid phase boundary) in our previous GCMC simulations [14] is more suitable for combination with a metadynamics study, since only one homogeneous phase exists in the confined region.
GCMC simulations have been considered a natural way and a standard method to simulate an open system with constant temperature T, volume V, and chemical potential µ (µVT system) imposed by an external reservoir while the number of particles is allowed to vary. Typical examples are clay swelling [28][29][30][31] and the quasistatic friction in aqueous film in an SFA experiment [32]. In GCMC simulation, at a given gap width between two confining walls, a hysteresis loop with two distinct branches (adsorption and desorption) is often observed and widely documented, which was not considered in our previous study [14]. How to determine the equilibrium transition within the hysteresis loop is a critical challenge. In this study, a recently proposed method called the mid-density scheme was applied to determine the equilibrium state in a hysteresis loop [33]. Compared with conventional but computationally expensive methods such as thermodynamic integration [34,35] and gauge cell method [36,37], the mid-density scheme is a simplified but useful approach to identify the equilibrium phase transition in a hysteresis loop.
The aim of this paper is to obtain an appropriate value of chemical potential extracted from hybrid grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations. The criterion is that it can correctly reproduce the thermodynamic states observed in LVMD simulation, rather than accurately predict the energetically favorable states. In other words, LVMD simulations can serve as a guide in searching the proper value of chemical potential adopted for a GCMC simulation. The equilibrium configuration obtained by this approach can be viewed as a reference and further used as the initial configuration in metadynamics to investigate the free energy profile relative to an order parameter at a given gap distance between confining solid walls.
The rest of the paper is organized as follows. In Section 2, we briefly describe the simulation models and the GCMC simulation techniques. Results are given in Section 3, followed by a summary and discussions in Section 4.

Simulation Models and Methods
As illustrated in Figure 1, argon atoms, as a nonpolar model lubricant, are confined between two rigid face-centered cubic (FCC) solid walls. A simple and computationally cheap Lennard-Jones (LJ) atomic pair potential is used to model the interactions between argon molecules according to the literature, i.e., σ = 0.3405 nm and ε = 0.2381 kcal/mol [38]. Following our previous publications [12,16], two types of contacts are considered: commensurate (Figure 1a,b) and incommensurate (Figure 1c,d). The two FCC crystal walls have the z direction along the [111] direction, while the x and y directions are along the 110 and 112 directions, respectively. For convenience, the relevant parameters are summarized in Table 1. For both contacts, the wall-fluid interaction between argon and solid wall is the same as the argon-argon interaction (i.e., ε wf = ε ff = ε, where w and f stand for "wall" and "fluid," respectively).
The rest of the paper is organized as follows. In Sections 2, we briefly describe the simulation models and the GCMC simulation techniques. Results are given in Sections 3, followed by a summary and discussions in Sections 4.

Simulation Models and Methods
As illustrated in Figure 1, argon atoms, as a nonpolar model lubricant, are confined between two rigid face-centered cubic (FCC) solid walls. A simple and computationally cheap Lennard-Jones (LJ) atomic pair potential is used to model the interactions between argon molecules according to the literature, i.e., σ = 0.3405 nm and ε = 0.2381 kcal/mol [38]. Following our previous publications [12,16], two types of contacts are considered: commensurate (Figure 1a,b) and incommensurate (Figure 1c,d). The two FCC crystal walls have the z direction along the [111] direction, while the x and y directions are along the 11 0 ·and 112 ·directions, respectively. For convenience, the relevant parameters are summarized in Table 1. For both contacts, the wall-fluid interaction between argon and solid wall is the same as the argon-argon interaction (i.e., εwf = εff = ε, where w and f stand for "wall" and "fluid," respectively).    Periodic boundary conditions (PBCs) are implemented in the x and y directions to approximate the infinite extent of lateral surfaces. All the simulations are carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [39]. A standard cutoff distance of 8.5 Å (about 2.5 σ) is used for the LJ interactions, and the time  [40,41] at a temperature of 85 K, corresponding to a liquid state of argon which is between the melting point (83.8 K) and the boiling point (87.3 K) of argon.
The local structural ordering of atoms is characterized by the OVITO visualization package [42] with the polyhedral template matching (PTM) method [43]. The PTM method can reliably identify the most common ordered crystalline structures of condensed phases ranging from face-centered cubic (FCC), hexagonal close-packed (HCP), body-centered cubic (BCC), simple cubic (SC), and icosahedral (ICO) to others (unknown coordination structure). It can still work effectively in the presence of strong thermal displacements and strains.
The GCMC simulation conducts exchanges of atoms or molecules with a virtual ideal gas reservoir. Thus, it is considered a natural way to simulate an open system characterized by fixed temperature T, volume V, and chemical potential µ. If all the inserted atoms in the confined region can also be moved using a regular NVT ensemble control (the volume and temperature are held constant), simulations in the grand canonical ensemble (µVT, constant chemical potential, constant volume, and constant temperature) can be performed. In such a hybrid GCMC/MD simulation, the inserted molecules can be well equilibrated.
Our hybrid GCMC/MD simulation attempts 1000 insertion/deletion moves and 1000 translational moves every 100 timesteps. In each simulation, at least 5 million timesteps are performed. The temperature of the imaginary reservoir is set to be equivalent to the target temperature used in NVT (T = 85 K) so that the imaginary reservoir is in thermal equilibrium with the simulation cell. When an argon atom is to be inserted, its coordinates are chosen at a random place within the confined region, and new atom velocities are randomly chosen from the specified temperature distribution given by T.
Hybrid GCMC/MD simulations can be initiated with an empty box or a box full of argon molecules, which will give the adsorption and desorption branches respectively. For the former, we continue to increase the chemical potential each time and the resulting configuration of the previous simulation is used as the starting configuration for the subsequent simulation. For the latter, we repeat the same procedure with decreasing chemical potential. In many microporous systems, the adsorption and desorption branches are distinct, exhibiting a hysteresis loop, with the equilibrium phase transition undetermined.
Compared with conventional but computationally expensive methods, the middensity scheme is a simplified but useful approach to identify the equilibrium phase transition in a hysteresis loop [33]. It has been successfully applied to study argon adsorption in graphite slit pores and in silica cylindrical pores [33,44]. The researchers who proposed this method also demonstrated [45] that its validity for the determination of the equilibrium state is supported by the commonly used canonical ensemble kinetic Monte Carlo (CE-kMC) [46,47] and gauge cell methods [36,37,48].
The detailed procedures are as follows: (1) For each chemical potential, there are two states: low-density state (with average N low argon molecules in the confined region) on the adsorption branch and high-density states (with average N high argon molecules in the confined region) on the desorption branch. (2) The number of particles in the mid-density state is the average of the number of particles in the low-density state and the high-density state, i.e., N mid = (N high + N low )/2. (3) The mid-density state can be constructed by randomly deleting (N high − N low )/2 molecules from the high-density state. (4) Run the NVT MD simulation to completely relax the system with N mid particles. (5) As a final step, the hybrid GCMC/MD simulation is performed for a sufficiently long time and the final state obtained is considered to be the stable state. It is summarized in the flowchart shown in Figure 2.

Results and Discussion
Our previous LVMD simulations show that the critical layer number, nc, at which the confinement-induced liquid-like to solid-like phase transition occurs, does not depend on the size of the contact area, but depends on the nature of contact [16]. As to the commensurate contact, this number is nc = 7, while in the case of the incommensurate contact, it is nc = 5, the same critical layer number for mica-cyclohexane-mica contact [15]. Our special interest is in the equilibrium state in both commensurate and incommensurate contacts at nc and nc + 1 monolayer thickness.

Commensurate Contact Simulations
For the perfect commensurate contact, hybrid GCMC/MD simulations are performed at n = 7 and n = 8 monolayer thickness. The adsorption-desorption isotherms are shown in Figure 3. The black and red curves represent adsorption and desorption branches respectively, which form two very broad hysteresis loops. The overall pattern of each adsorption and desorption branch is a three-step shape, representing the gas, liquid and solid phases respectively.

Results and Discussion
Our previous LVMD simulations show that the critical layer number, n c , at which the confinement-induced liquid-like to solid-like phase transition occurs, does not depend on the size of the contact area, but depends on the nature of contact [16]. As to the commensurate contact, this number is n c = 7, while in the case of the incommensurate contact, it is n c = 5, the same critical layer number for mica-cyclohexane-mica contact [15]. Our special interest is in the equilibrium state in both commensurate and incommensurate contacts at n c and n c + 1 monolayer thickness.

Commensurate Contact Simulations
For the perfect commensurate contact, hybrid GCMC/MD simulations are performed at n = 7 and n = 8 monolayer thickness. The adsorption-desorption isotherms are shown in Figure 3. The black and red curves represent adsorption and desorption branches respectively, which form two very broad hysteresis loops. The overall pattern of each adsorption and desorption branch is a three-step shape, representing the gas, liquid and solid phases respectively.
As to n = 7 monolayer, the adsorption branch starts with µ = −2.40 kcal/mol, in which the system is in a state of dilute gas phase and the argon atoms are moving close to the confining walls ( Figure 4a). As the chemical potential increases, the number of atoms in the gas phase increases and two contact layers are formed (Figure 4b). At µ = −2.225 kcal/mol, the system experiences a gas to liquid phase transition with a sharp increase in the number of particles (Figure 4c). At µ = −2.19 kcal/mol, the system undergoes a liquid to solid phase transition, in which FCC layers and HCP layers alternate to form 7 solid layers ( Figure 4d). The desorption branch starts with µ = −2.15 kcal/mol, in which the solidified film is formed in a polycrystalline state that has FCC and HCP crystal structures. The grain boundaries between FCC and HCP crystals are clearly seen in the two outmost contact layers as shown in Figure 4e. The middle three layers are all in HCP structures. Due to existence of grain boundaries and defects, the number of argon particles in such a polycrystalline state is smaller than that in solid state in the adsorption branch in which the packing structure is more compact (Figure 4d). As the chemical potential decreases, the number of atoms in solid state remains constant and the microstructure does not change a lot. At µ = −2.325 kcal/mol, the system experiences a solid (Figure 4f) to liquid phase transition (Figure 4g). At µ = −2.40 kcal/mol, the system goes back to the gas phase (Figure 4h), completing the hysteresis loop.
surate contact, this number is nc = 7, while in the case of the incommensurate conta nc = 5, the same critical layer number for mica-cyclohexane-mica contact [15]. Our s interest is in the equilibrium state in both commensurate and incommensurate cont nc and nc + 1 monolayer thickness.

Commensurate Contact Simulations
For the perfect commensurate contact, hybrid GCMC/MD simulations are perfo at n = 7 and n = 8 monolayer thickness. The adsorption-desorption isotherms are s in Figure 3. The black and red curves represent adsorption and desorption branch spectively, which form two very broad hysteresis loops. The overall pattern of ea sorption and desorption branch is a three-step shape, representing the gas, liqui solid phases respectively.   As to n = 7 monolayer, the adsorption branch starts with μ = −2.40 kcal/mol, in which the system is in a state of dilute gas phase and the argon atoms are moving close to the confining walls ( Figure 4a). As the chemical potential increases, the number of atoms in the gas phase increases and two contact layers are formed (Figure 4b). At μ = −2.225 kcal/mol, the system experiences a gas to liquid phase transition with a sharp increase in the number of particles (Figure 4c). At μ = −2.19 kcal/mol, the system undergoes a liquid to solid phase transition, in which FCC layers and HCP layers alternate to form 7 solid layers ( Figure 4d). The desorption branch starts with μ = −2.15 kcal/mol, in which the solidified film is formed in a polycrystalline state that has FCC and HCP crystal structures. The grain boundaries between FCC and HCP crystals are clearly seen in the two outmost contact layers as shown in Figure 4e. The middle three layers are all in HCP structures. Due to existence of grain boundaries and defects, the number of argon particles in such a polycrystalline state is smaller than that in solid state in the adsorption branch in which the packing structure is more compact (Figure 4d). As the chemical potential decreases, the number of atoms in solid state remains constant and the microstructure does not change a lot. At μ = −2.325 kcal/mol, the system experiences a solid (Figure 4f) to liquid phase transition (Figure 4g). At μ = −2.40 kcal/mol, the system goes back to the gas phase (Figure 4h), completing the hysteresis loop. As the next step, the mid-scheme is applied to these adsorption and desorption isotherms to obtain the equilibrium phase as shown in the blue curves in Figure 3a. When μ ≤ −2.30 kcal/mol, the system is in the gas phase and the number particles in the gas phase is notably larger than that in the adsorption branch. When chemical potential is in the range of [−2.29, −2.22] kcal/mol, the equilibrium state is the liquid phase. As can be seen, the solid state exists when chemical potential is larger than −2.225 kcal/mol. Still, two possible solid states are observed: a polycrystalline state similar to that in Figure 4e with smaller chemical potential and 7 solid layers with stacking fault similar to that in Figure  4d corresponding to larger chemical potential.
We repeat the same procedure for the case of n = 8 monolayer thickness. The shape of the hysteresis loop is similar to the n = 7 monolayer case and, as expected, it is wider. In this case, there is no much difference between the two solid states in adsorption and desorption branches (Figure 5d,e). They are both in the polycrystalline state and the grain As the next step, the mid-scheme is applied to these adsorption and desorption isotherms to obtain the equilibrium phase as shown in the blue curves in Figure 3a. When µ ≤ −2.30 kcal/mol, the system is in the gas phase and the number particles in the gas phase is notably larger than that in the adsorption branch. When chemical potential is in the range of [−2.29, −2.22] kcal/mol, the equilibrium state is the liquid phase. As can be seen, the solid state exists when chemical potential is larger than −2.225 kcal/mol. Still, two possible solid states are observed: a polycrystalline state similar to that in Figure 4e with smaller chemical potential and 7 solid layers with stacking fault similar to that in Figure 4d corresponding to larger chemical potential. We repeat the same procedure for the case of n = 8 monolayer thickness. The shape of the hysteresis loop is similar to the n = 7 monolayer case and, as expected, it is wider. In this case, there is no much difference between the two solid states in adsorption and desorption branches (Figure 5d,e). They are both in the polycrystalline state and the grain boundaries between FCC and HCP crystals are in the two or three outmost contact layers. Based on the mid-density scheme applied to the isotherms, the liquid state is the stable phase when chemical potential is smaller than −2.20 kcal/mol but larger than −2.29 kcal/mol. boundaries between FCC and HCP crystals are in the two or three outmost contact layers. Based on the mid-density scheme applied to the isotherms, the liquid state is the stable phase when chemical potential is smaller than −2.20 kcal/mol but larger than −2.29 kcal/mol. Since the critical layer number for the perfect commensurate contact is nc = 7 in LVMD simulations, the confined film is in liquid state at n = 8 layers and is in solid state at n = 7 layers. Taken together, the appropriate chemical potential should be in the range of [−2.225, −2.20] kcal/mol, which can reproduce the correct equilibrium phase for n = 7 and 8 layers observed in LVMD simulations.
In the next step, we select chemical potential μ = −2.21 kcal/mol for the rest of simulations with varying thickness. As can be seen in Table 2, all the equilibrium phases predicted in GCMC simulations ( Figure 6) are consistent with the observations of LVMD simulations (see Figure A1), which confirms that μ = −2.21 kcal/mol is an appropriate choice for running GCMC/MD simulations. With this chemical potential, GCMC/MD can produce equilibrium configuration for future simulations coupled with metadynamics. Table 2. The equilibrium phase in different thicknesses predicted by GCMC/MD simulations with chemical potential μ = −2.21 kcal/mol ( Figure 6) and liquid-vapor molecular dynamics (LVMD) simulations ( Figure A1).  Since the critical layer number for the perfect commensurate contact is n c = 7 in LVMD simulations, the confined film is in liquid state at n = 8 layers and is in solid state at n = 7 layers. Taken together, the appropriate chemical potential should be in the range of [−2.225, −2.20] kcal/mol, which can reproduce the correct equilibrium phase for n = 7 and 8 layers observed in LVMD simulations.
In the next step, we select chemical potential µ = −2.21 kcal/mol for the rest of simulations with varying thickness. As can be seen in Table 2, all the equilibrium phases predicted in GCMC simulations ( Figure 6) are consistent with the observations of LVMD simulations (see Figure A1), which confirms that µ = −2.21 kcal/mol is an appropriate choice for running GCMC/MD simulations. With this chemical potential, GCMC/MD can produce equilibrium configuration for future simulations coupled with metadynamics. Table 2. The equilibrium phase in different thicknesses predicted by GCMC/MD simulations with chemical potential µ = −2.21 kcal/mol ( Figure 6) and liquid-vapor molecular dynamics (LVMD) simulations ( Figure A1).

Incommensurate Contact Simulations
We now consider the incommensurate contact between the confining wall and the solidified film. The critical number of monolayers in the film for liquid-like to solid-like phase transition is nc = 5. Following the same procedure in the case of commensurate contact, the hybrid GCMC/MD simulations are performed at n = 5 and n = 6 monolayer thickness. The adsorption and desorption isotherms are shown in Figure 7. The black and red curves represent adsorption and the desorption branches, respectively. The blue curves show the equilibrium states by applying the mid-density scheme to the adsorption-desorption isotherms. The overall pattern of the hysteresis loop is similar to that in the commensurate case except the desorption branch at n = 5 monolayer in which the liquid phase is absent.

Incommensurate Contact Simulations
We now consider the incommensurate contact between the confining wall and the solidified film. The critical number of monolayers in the film for liquid-like to solid-like phase transition is n c = 5. Following the same procedure in the case of commensurate contact, the hybrid GCMC/MD simulations are performed at n = 5 and n = 6 monolayer thickness. The adsorption and desorption isotherms are shown in Figure 7. The black and red curves represent adsorption and the desorption branches, respectively. The blue curves show the equilibrium states by applying the mid-density scheme to the adsorptiondesorption isotherms. The overall pattern of the hysteresis loop is similar to that in the commensurate case except the desorption branch at n = 5 monolayer in which the liquid phase is absent.

Incommensurate Contact Simulations
We now consider the incommensurate contact between the confining wall and solidified film. The critical number of monolayers in the film for liquid-like to solid phase transition is nc = 5. Following the same procedure in the case of commensurate tact, the hybrid GCMC/MD simulations are performed at n = 5 and n = 6 monolayer th ness. The adsorption and desorption isotherms are shown in Figure 7. The black and curves represent adsorption and the desorption branches, respectively. The blue cu show the equilibrium states by applying the mid-density scheme to the adsorption sorption isotherms. The overall pattern of the hysteresis loop is similar to that in the c mensurate case except the desorption branch at n = 5 monolayer in which the liquid p is absent.   As for the n = 5 monolayer, in the adsorption branch, the equilibrium state is the gas phase with chemical potential in the range of [−2.60, −2.33] kcal/mol (Figure 8a). Liquid state (Figure 8b) can only exist in a very narrow range (from −2.325 kcal/mol to −2.320 kcal/mol), which is followed by solid phase when chemical potential µ ≥ −2.31. It should be noted that the normal of the (111) close packed plane of the solid structure is nearly along the diagonal direction in the lateral x-y plane (Figure 8c). The desorption branch starts with µ = −2.15 kcal/mol, in which the solidified film is formed. Compared with the solidified structure in the adsorption branch (Figure 8c), the lattice orientation is different and its close packed normal is off the diagonal of the lateral x-y plane (Figure 8d). The slightly different configurations can explain the difference in the number of particles in solidified structures (Figure 8c,d). As the chemical potential decreases, the number of atoms in solid state remains constant and the microstructure does not change a lot (Figure 8e). At µ = −2.515 kcal/mol, the system experiences a solid (Figure 8e) to gas phase transition (Figure 8f), completing a hysteresis loop. In this case, no liquid phase is observed. In the next step, the mid-scheme is applied to these adsorption and desorption isotherms to obtain the equilibrium phase as shown in the blue curves in Figure 7a. It can be seen that the equilibrium phase is the solid state when chemical potential is larger than −2.31 kcal/mol. For the case of n = 6 monolayer, even though the lattice orientations of the solidified structures in the adsorption and desorption branches are different (Figure 9c,d), there is no distinct difference in the number of particles in the solid structures. As can be seen, the liquid state is the stable phase when chemical potential is in the range of [−2.33, −2.25] kcal/mol. Considering the critical layer number for the incommensurate contact being n c = 5 in LVMD simulations, the appropriate chemical potential should be in the range of [−2.31, −2.25] kcal/mol, which can reproduce the correct equilibrium phase for n = 5 (in solid state) and n = 6 (in liquid state) monolayer observed in LVMD simulations. In the following, chemical potential µ = −2.28 kcal/mol is adopted for the rest of the simulations with varying thicknesses as shown in Table 3. All the equilibrium phases predicted in GCMC simulations are consistent with the observations of LVMD simulations, confirming that µ = −2.28 kcal/mol is an appropriate choice.
Coatings 2021, 11, x FOR PEER REVIEW 9 of 14 As for the n = 5 monolayer, in the adsorption branch, the equilibrium state is the gas phase with chemical potential in the range of [−2.60, −2.33] kcal/mol (Figure 8a). Liquid state (Figure 8b) can only exist in a very narrow range (from −2.325 kcal/mol to −2.320 kcal/mol), which is followed by solid phase when chemical potential μ ≥ −2.31. It should be noted that the normal of the (111) close packed plane of the solid structure is nearly along the diagonal direction in the lateral x-y plane (Figure 8c). The desorption branch starts with μ = −2.15 kcal/mol, in which the solidified film is formed. Compared with the solidified structure in the adsorption branch (Figure 8c), the lattice orientation is different and its close packed normal is off the diagonal of the lateral x-y plane (Figure 8d). The slightly different configurations can explain the difference in the number of particles in solidified structures (Figure 8c,d). As the chemical potential decreases, the number of atoms in solid state remains constant and the microstructure does not change a lot (Figure  8e). At μ = −2.515 kcal/mol, the system experiences a solid (Figure 8e) to gas phase transition (Figure 8f), completing a hysteresis loop. In this case, no liquid phase is observed. In the next step, the mid-scheme is applied to these adsorption and desorption isotherms to obtain the equilibrium phase as shown in the blue curves in Figure 7a. It can be seen that the equilibrium phase is the solid state when chemical potential is larger than −2.31 kcal/mol. For the case of n = 6 monolayer, even though the lattice orientations of the solidified structures in the adsorption and desorption branches are different (Figure 9c,d), there is no distinct difference in the number of particles in the solid structures. As can be seen, the liquid state is the stable phase when chemical potential is in the range of [−2.33, −2.25] kcal/mol. Considering the critical layer number for the incommensurate contact being nc = 5 in LVMD simulations, the appropriate chemical potential should be in the range of [−2.31, −2.25] kcal/mol, which can reproduce the correct equilibrium phase for n = 5 (in solid state) and n = 6 (in liquid state) monolayer observed in LVMD simulations. In the following, chemical potential μ = −2.28 kcal/mol is adopted for the rest of the simulations with varying thicknesses as shown in Table 3. All the equilibrium phases predicted in GCMC simulations are consistent with the observations of LVMD simulations, confirming that μ = −2.28 kcal/mol is an appropriate choice.   Figure A2.  Figure A2.

Summary and Further Discussion
It has been discussed that LVMD ensemble is not suitable for metadynamics simulations, because of the coexistence of multi-phases. Instead, simple geometry with slit pores can be a proper simulation setup. In this work, hybrid GCMC/MD simulations are carried out to investigate the adsorption and desorption isotherms of argon molecules confined between commensurate and incommensurate contacts in nanoscale thickness. By applying the mid-density scheme to the hysteresis loops, the equilibrium structures associated with n c and n c + 1 layers can be obtained. The appropriate chemical potentials can be determined if the equilibrium structures predicted by GCMC/MD simulations are consistent with those observed in LVMD simulations. Such equilibrium structures at varying thicknesses obtained using GCMC/MD simulations can be used as reasonable initial configurations for further metadynamics studies to investigate the free energy profile as a function of structural order parameters.
The shift of the chemical potential due to the presence of nanoconfinement makes the selection of chemical potential µ difficult, and in many situations the choice is merely based on the bulk chemical potential. In our previous study [14], we discussed that the bulk chemical potential µ for argon at around 85 K obtained either from Widom's insertion method in MC calculations or from the integral equation theory [49] cannot be used in GCMC simulation. In this study, we find that the appropriate chemical potential is in the range of [−2.225, 2.20] kcal/mol for the commensurate contact and [−2.31, −2.25] kcal/mol for the incommensurate contact. There is even no overlap between them, which means that no unique chemical potential can be adopted for both cases. Given the fact that confinement usually induces a shift in chemical potential when considering the equilibrium between the slit pore and bulk fluid, we expect that the shift due to the different nature of contact should also be different.
It should be noted that in previous works, researchers mainly focus on the capillary condensation (the gas-liquid transition) in the pore, in which only two phases, gas and liquid, are involved in the GCMC simulations. When the confined space comes down to the boundary lubrication (BL) regime, three phases come into play, which makes the hysteresis loop more complicated and compounds the problems of determining the equilibrium state. As a first attempt, we apply the mid-density scheme to such a hysteresis loop. There is still some uncertainty about whether the mid-density scheme can be applied to the solid phase. However, our purpose is not to construct equilibrium phase transition in a very accurate and rigorous manner, but to construct a reasonable initial configuration at a given gap distance by GCMC/MD simulations to further investigate the free energy profile vs. structural order parameters. To accurately identify the location of the phase equilibrium, the gauge cell MC method [36,37] needs to be applied in future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Equilibrium Phase in LVMD Simulations
The equilibrium LVMD simulations of argon molecules confined between commensurate and incommensurate contacts are performed at a certain wall distance. It is found that the equilibrium configurations observed in LVMD simulations are consistent with those observed in GCMC/MD simulations. It should be noted that in this work, the computational cost of both LVMD and GCMC/MD simulations are within the same order of magnitude (every single LVMD or GCMC/MD simulation at each gap thickness is about a few hours long).
Coatings 2021, 11, x FOR PEER REVIEW 12 of 14 that the equilibrium configurations observed in LVMD simulations are consistent with those observed in GCMC/MD simulations. It should be noted that in this work, the computational cost of both LVMD and GCMC/MD simulations are within the same order of magnitude (every single LVMD or GCMC/MD simulation at each gap thickness is about a few hours long).       Figure 10.