THE JOURNAL OF CHEMICAL PHYSICS 142, 074701 (2015)

Improved modeling of two-dimensional transitions in dense phases on crystalline surfaces. Krypton–graphite system E. A. Ustinova) Ioffe Institute, 26 Polytechnicheskaya, St. Petersburg 194021, Russia

(Received 22 December 2014; accepted 2 February 2015; published online 17 February 2015) This paper presents a refined technique to describe two-dimensional phase transitions in dense fluids adsorbed on a crystalline surface. Prediction of parameters of 2D liquid–solid equilibrium is known to be an extremely challenging problem, which is mainly due to a small difference in thermodynamic functions of coexisting phases and lack of accuracy of numerical experiments in case of their high density. This is a serious limitation of various attempts to circumvent this problem. To improve this situation, a new methodology based on the kinetic Monte Carlo method was applied. The methodology involves analysis of equilibrium gas–liquid and gas–solid systems undergoing an external potential, which allows gradual shifting parameters of the phase coexistence. The interrelation of the chemical potential and tangential pressure for each system is then treated with the Gibbs–Duhem equation to obtain the point of intersection corresponding to the liquid/solid–solid equilibrium coexistence. The methodology is demonstrated on the krypton–graphite system below and above the 2D critical temperature. Using experimental data on the liquid–solid and the commensurate–incommensurate transitions in the krypton monolayer derived from adsorption isotherms, the Kr–graphite Lennard–Jones parameters have been corrected resulting in a higher periodic potential modulation. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4908035]

I. INTRODUCTION

Liquid–solid transition attracts much attention in the literature for several decades. In part, this is motivated by the difficulty of evaluating the true thermodynamic transition, which is a consequence of relatively small difference in the density and thermodynamic potentials of coexisting phases. For this reason, so-called “direct” approaches, which explicitly involve the solid–liquid interface1 or use gradual change in temperature or pressure to detect the abrupt change in density or structure of the simulated phase,2–4 usually encounter the problem of its superheating or supercooling and, therefore, do not lead to rigorous estimation of the transition temperature. Such approaches as the hysteresis method4 or voids method3 are based on some empirical assumptions, which do not provide a reliable estimation of the equilibrium melting/freezing transition. To overcome this difficulty, various techniques based on thermodynamic integration have been suggested. Thus, one can mention the so-called reference state path.5,6 The method involves an artificial solid, usually the Einstein crystal, using an external potential to ensure a reversible phase transition, which allows to calculate the free energy of the real crystal. The reference ideal gas is taken to determine the free energy of liquid. The equivalence of the specific Gibbs free energy for the solid and liquid corresponds to the liquid–solid thermodynamic equilibrium. Once a point of the phase coexistence is known, a)Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2015/142(7)/074701/10/$30.00

the melting curve can be further evaluated by integration of the Gibbs–Duhem equation.7,8 In other group of methods, a similar idea is used to determine directly the difference between free energies of the liquid and solid without invoking the reference system.9 In both cases, the λ-integration and Gaussian wells are used as a constraint to artificially keep the fluid in the crystalline form over intermediate stages. The review of these approaches in case of 3D systems can be found elsewhere.10–13 Although methods based on the free energy difference evaluation are generally more rigorous than the “direct” approaches, one can mention the recently published “interface pinning” method,14,15 which seems to overcome the problem of instability of the liquid–solid system by introducing a harmonic field between the two phases. Two-dimensional phase transitions have some distinctions relative to their bulk counterpart mainly due to the confinement effect in micropores and the effect of external potential exerted by an open surface that stabilizes the 2D layer. There are two different situations. In the first case, when the surface is not completely filled by the adsorbed fluid, the solid or liquid molecular layer formed at a sufficiently low temperature is in equilibrium with the gas phase. Therefore, the liquid–solid transition occurs at the tangential pressure close to zero. If the surface is completely filled, the tangential pressure can achieve several thousand bars, which significantly increases the melting/freezing temperature. For this reason, it is convenient to consider these two cases separately. One of the most interesting and relatively simple systems is krypton adsorbed on graphite where both the liquid–solid and solid–solid (commensurate–incommensurate

142, 074701-1

© 2015 AIP Publishing LLC

074701-2

E. A. Ustinov

(C-IC)) transitions can be observed. This system has been intensively investigated experimentally16–24 and theoretically using molecular simulations.25–30 All methods developed for bulk phases can be applied to analysis of 2D phase transitions. Thus, Das and Singh31 used the pseudo-supercritical transformation path of Grochola9 to determine the liquid–solid phase transition in slit micropores. One of main difficulties in analysis of phase equilibrium involving the crystalline phase is uncertainty of its lattice constant at a given chemical potential (or pressure) and temperature. Simulation in the box of constant arbitrarily chosen volume models an artificial system which can be thermodynamically self-consistent but not corresponding to the real crystalline phase at the same pressure and temperature. This problem arises from disappearance of the explicit link between the chemical potential and the pressure in the crystal because at a given temperature, the latter is nearly the singlevalued function of the lattice constant over a wide range of chemical potential. This issue has been addressed by Sweatman et al.32 who noted that even the Gibbs ensemble does not solve the problem. This difficulty was shown to circumvent by adjusting the simulation box volume in the course of gradual increase of the chemical potential so that the corresponding change in the tangential pressure obeys the Gibbs–Duhem equation.33,34 In order to determine parameters of the 2D liquid–solid transition, we have developed a method34,35 similar to the Maxwell rule of equal areas, in which the van der Waals loop was presented as a dependence of the surface density on the tangential pressure. This method proved to be highly reliable in the case of structureless substrate. However, in the case of the Kr–graphite system, the method does not work because the tangential pressure cannot be determined in the uniform 2D crystal commensurate with the graphite. This is because the total tangential pressure is shared by the crystal and graphite surface, with the commensurate krypton lattice being stretched resulting in the graphite compression. The latter occurs only in the krypton gas/liquid–solid interface and, therefore, can be determined if the interface is explicitly involved into consideration. It has been done in the framework of a new methodology,36 where the 2D gas–solid (commensurate or incommensurate) system on the crystalline substrate is simulated in the same elongated box. There are three advantages of such a technique: (i) the tangential pressure in x y-plane parallel to the surface which acts on the substrate can be determined by summation of resulting forces between carbon atoms and adsorbed molecules in the transition zone (outside this zone, the resulting forces are zero on average); (ii) the chemical potential of the whole equilibrium system (and, consequently, that of the crystalline phase) is determined accurately due to existence of the gas phase and the transition zone; (iii) the average distance between neighboring molecules in the ordered phase along the normal to the interface is not constrained by the simulation box making the lattice constant a self-adjusting parameter. The simulation box dimension parallel to the gas–solid interface is adjusted to equilibrate components of the tangential pressure in the crystalline phase. Additionally, an external potential is imposed on the gas phase of the system to shift the equilibrium toward higher or lower values of the chemical potential and, hence, the tangential

J. Chem. Phys. 142, 074701 (2015)

pressure which allows one to find the point where both variables of different dense phases are equal.36 The methodology proved to be quite successful as applied to the C-IC transition. The aim of this paper is to extent the methodology to thermodynamic analysis of the 2D liquid–solid transition at temperatures below and above the 2D krypton critical point and to correct the Lennard–Jones (LJ) parameters of carbon atoms using experimental data on the liquid–solid transition24 and commensurate–incommensurate transition.37 This is a very important point because the widely used semi-empirical Steele parameters for carbon atom38 (i.e., 0.34 nm and 28 K for the collision diameter, σss, and the potential, ε ss/k B, respectively) still do not have a reliable experimental confirmation, though there is some evidence that the collision diameter of carbon atom is lower than 0.34 nm providing a higher periodic modulation of the potential exerted by the graphite surface.39,40 This also means that kernels used in the pore size distribution analysis of nanoporous materials by argon and nitrogen adsorption isotherms should be recalculated. II. METHODOLOGY AND SIMULATION DETAILS

The kinetic Monte Carlo method41 (kMC) was used in simulation of the 2D gas–solid/liquid system in a canonical ensemble. A feature of the kMC scheme is that if the displaced molecule A overlaps with the molecule B, both molecules take very high positive potential energies and, therefore, high mobility (rates). Then, according to the kMC scheme, the two molecules have a comparable probability to be chosen for the next MC step, while in the standard Metropolis algorithm, this displacement is not accepted. Therefore, in the case of dense phases, kMC provides much faster change of configurations. An additional advantage of the kMC as applied to the gas–solid/liquid coexistence is that the gas phase thermodynamic properties are determined with extremely high reliability owing to a large mobility of gas molecules providing a good statistics despite the low density of this phase. In its turn, this permits a reliable evaluation of the chemical potential of the coexisting dense phase. At the same time, it should be mentioned a small deficiency of the kMC appearing in systems which include the coexisting crystalline and gas phases in the same simulation box. In this case, the most intensive change occurs in the transition zone which combines a high mobility of molecules and relatively high density. This is a reason of relatively small number of molecular displacements in the crystalline phase. In order to balance statistics in coexisting phases and the transition zone, the standard Metropolis algorithm has been incorporated into the kMC scheme. Such a combined scheme implies that in each 1000 kMC step, a series of 1000 Metropolis MC steps is accomplished, with the maximum displacement being adjusted to keep the acceptance ratio close to 0.5. This stage has been included to the kMC scheme only to provide an appropriate small rearrangement of the crystalline phase to speed up a necessary change of its lattice constant, but it is omitted in the averaging procedure, so that the chemical potential, pressure, and internal energy are still determined solely in the framework of kMC.

074701-3

E. A. Ustinov

Simulations were performed in elongated rectangle box with the dimensions ratio L x /L y = 4(31/2/2) and L z = 5σff (normal to the surface). The 12-6 LJ parameters σff and ε ff /k B for Kr–Kr interaction were taken as 0.36 nm and 171 K. The dense phase was kept at the center of the box and occupied one-half of its area in x y-plane parallel to the graphite. Eight hundred molecules were taken for simulation of the gas–solid equilibrium. In the case of gas–liquid coexistence, the number of molecules in each run was chosen during the equilibration stage so that the liquid phase occupied approximately onehalf of the x y-plane. Simulation of the gas–commensurate solid and the gas–liquid equilibrium was carried out in the box having L y dimension of n(31/2/2)3a0, where n = 20 and a0 = 0.142 nm (the distance between neighboring carbon atoms in the graphite lattice). The symmetry axis of graphite was parallel to the L x dimension. If the crystalline phase was incommensurate, the value of L y dimension was adjusted during the equilibration stage to enforce equality of x and y projections of the tangential pressure. The periodic boundary conditions were imposed only along x and y directions. The number of kMC steps in all runs was 25 × 106 for equilibration and 50 × 106 for sampling. The contribution of the outermost graphite layer to the solid–fluid potential was modeled by a direct summation of the C–Kr LJ pair potentials within the circle of radius 10a0 from the projection of the krypton atom on the graphite surface, which involves about 140 carbon atoms.36 The contribution of all rest carbon atoms of the same layer was evaluated by integration of the 12-6 LJ potential over the x y plane for r > 10a0. The potential exerted by the second and other graphite layers (up to 10) was modeled as a sum of 10-4 potentials, with the spacing being equal to 0.335 nm. The LJ parameters for the krypton–carbon atom interaction can be estimated with the Lorentz–Berthelot mixing rule using the C–C parameters defined by Steele38 as 0.34 nm and 28 K. This gives 0.35 nm for the C–Kr collision diameter, σ s f , and 69.2 K for the potential, ε s f /k B. These parameters, however, can be adjusted based on experimental data on phase transitions in the Kr–graphite system.24,37 As is shown previously,36 the gas–liquid/solid equilibrium can be shifted towards higher or lower value of the chemical potential (and, hence, the tangential pressure) by imposing an external potential on the gas phase. A simple form of the potential based on a polynomial of the third order can be expressed as follows:36

J. Chem. Phys. 142, 074701 (2015)

be positive or negative. Therefore, variation of um allows one to evaluate the dependence of tangential pressure, pT , on the chemical potential, µ, in the condensed phase. The point of intersection of the pT − µ curves plotted for the gas–liquid and gas–solid systems corresponds to the liquid–solid coexistence. A. Tangential pressure

It is remarkable that the external potential offers a simple way to determine the difference between the tangential pressures in coexisting phases. The transition zone of the potential presents a potential barrier, so any molecule in that zone undergoes the force along x-direction which contributes to the overall pressure difference between the gas and the liquid (or solid) phases. The gas phase pressure is usually very small and can be determined quite reliably, so the sum of the latter and the pressure difference due to the potential barrier is the tangential pressure in the condense phase. The tangential pressure difference at a distance z from the surface is written as follows:   ext  1 1 du (x i ) . (3) ∆pext (z) = x dx 2 L x ∆zσ 3 ff

z

Improved modeling of two-dimensional transitions in dense phases on crystalline surfaces. Krypton-graphite system.

This paper presents a refined technique to describe two-dimensional phase transitions in dense fluids adsorbed on a crystalline surface. Prediction of...
2MB Sizes 0 Downloads 6 Views