Computer Simulation of Metastable Fluid States in the Lennard-Jones System

Using the Monte Carlo method in statistical mechanics, we have simulated high density metastable states. We find that nucleation from a three dimensional fluid array to a crystalline solid is possible, but that periodic boundary conditions and the small size of the system inhibit the formation of perfect crystals. Evidence for the existence of an amorphous solid state has also been observed, and the pair correlation function of this state exhibits some of the features associated with random close-packed arrays of hard spheres. The possible relation between these simulations and the formation of glassy states in real systems is briefly discussed.


Introduction
The rmal and s truc tural properties of me tastable fluid s tates associated with subcoolin g are of inte rest in several proble m s, including glass formation, nu cleation and c rystal growth , and di sordered alloys.
W e re port brie fly on Monte Carlo [1)1 simulations of metastable states in sys te ms of 108 and 256 particles interacting through the Lennard-Jones 6-12 potential function. The en e r gy and length parameters are denote d by E and u, r es pectively. With N A the Avogadro cons tant and k B the Boltzm a nn cons tant, the following r educed units are used: volume, u 3 ; press ure, E/u 3 ; energy, N AE; and te mper ature, El k B • W e r e port simulations at t emperatures of T= 0.80,1.17, and 2.74 for de nsitie s be tween 1.05 and 1.38. Press ures for these st ates are w ell into the kilo bar range . The computations proceed as follows. The coordinates of a di sorder e d fluid array, obtaine d from an equilibrate d simulation at a liquid-like de ns ity, are scale d to a high d e nsity at a te mper ature a t whic h the lowes t e ner gy state is that of a crystallin e (fcc) solid [2]. (The initial configuration thu s has the charac ter of a m etas table fluid state, whic h c an be realized in a I Figures in brack ets in dicate the Iit eratu_re refere nces al th e end of thi s paper. 59 re al fluid by unde rcoolin g or o ver co mpression. ) W e the n de termin e wh e the r the simula tion conver ges to partic ular valu es of th e press ure and internal e ner gy, and whe ther, on the molec ular le vel, it converges to a sta te whic h is not noticeably labile.
W e find that the therm odynamic prope rties do in fac t conver ge , afte r runs of the orde r of 10 7 confi gurations , but that the re sulting press ures and e nergies are higher than those of the crystalline solid at the same density and te mperature. Examination of the molecular struc ture indicates that some of these runs conve rged to dis torte d fcc crys tals , while others converged to comparatively amorphous s ta tes, with structures similar to the random close-pac ked s tructures studi ed by Finne y [3] and Bern al [4]. S uc h stru ctures have been disc ussed by Turnbull [5] a s a possible model for the ide al glass state.
Although the sequen ce of confi g uration s generated in a Monte Carlo simul ation is produ ced by a rando m process, each confi guration differ s from the precedin g one by a t mos t a s mall di s place m ent of a sin gle particle. H e nce it is not unre asona ble to consid er that th e progress of the simulation bears an a pproximate re lation to the progr ess of real tim e . Examin a tion of the molecular stru ctures r es ultin g fro m these simula tions may th er efore provid e so m e in sight into cr ystal nuclea-tion and glass formation in metastable fluids composed of simple molec ules. Nucleation phenomena similar to those reported h ere have recently been obse rved in molec ular dynamics simulations of the same model [6]. Previous work on soft-sphere systems [7] also indicates similar results for repulsive pote ntials whi ch vary as the inverse twelfth power of the internuclear separation.

Results and Discussion
The results for the press ure and en ergy at p =1.375 and T= 2.74 are shown in figure la; this density is -10 p 80 5 10 15 20 CONFIGURATIONS. 10 6 25 particles appeared to be locked into a fixed array undergoin g only small displacements from their mean positions as the co mputation proceeded. The press ure and energy converged to values of 80 and -0.5, respectively. The press ure and energy obtained from a si mulation of an fcc array at the same de nsity and temperature are 72 and -1.9, respec tively.
The pair correlation function, g(r) at 5 X 10 6 and at 27 X 10 6 is shown by the solid curves in figure   The das hed curves are ave rages ove r the most recent 105 configurations and the solid cu rv es is the cumul ative average. The breaks in the solid curve indicate th at the particle coordinates at th ese po int s we re used to initiate a nt: w simulation. about 30 percent beyond the freezing den sity whi ch was reported [2] for the T= 2.74 isotherm. The solid c urves ar e the c umulative averages and the dashed curves are the averages over the most recent 10 5 configurations. The breaks in the solid c urve at 9 X 10 6 and 16 X 10 6 configurations indicate that the particle coordinates at these points were used to initiate a new simulation. The very close agreement between the solid and dashed c urves beyond 17 X 10 6 clearly indicates convergence. (The corresponding results for the internal energy be have similarly.) The nearly horizontal curves at the left of fi gures la and 2a were characterized by exten sive displacements of the particles from their initial positions -in other words " liquid-like" behavior. As the pressure and energy decreased, the net particle displacements over each 10 5 configurations likewise decreased_ At the horizontal portion of the c urves at the end of each run, the molecular structure was no longer labile; that is, the and that the initially metastable fluid has converge d to an imperfect crystalline array_ Visual inspection of three-dimensional structure at 27 X 10 6 configurations by means of stereoscopic pairs of 35 mm photographs confirms the exis ten ce of distorted crystalline planes_ More striking evidence for nucleation to a crystalline array was obtained at p=1.105 and T=0.80, whi c h is less than one-h alf the fre ezing temperature [2] at this density. The results for the pressure and the average root mean square displace me nt of a particle (rms) from its initial position are shown in figure 2a.

Th e rm s (o pen circl es) is in lI n il s of Ih e disl ance parameter, u.
The press ure over the most recent 10 5 configuration converged to 8.93 and th e res ult for a run begun from fcc positions was 9.03. The intern al en ergy for both run s converged to th e same valu e, -7.17. Th e be havi or of the rm s s uggests th a t as the syste m equilibrates th e positions of th e atom s are not partic ularly labile. This can be qualitatively inte rpre ted as restric ted diffu sion. The pair correlation for the run s are shown in fi gure 2b , as before the das hed c urve is for the run begun 2.0 X 10 6 confi gurati ons in th e course of the run begun from the random array . Evide ntly th e local order in th e ra ndo m sys te m underwent dis tin ct ch anges be twee n 2.0 and 20 X 10 6 confi gura ti ons. The form ation of a pe rfec t c rys tallin e lattice in th ese simulations is greatly hind e red by the pe ri odi c boundary condition s whi c h ar e used to avoid surface e ffects. Unless the crystal nucleates with its prin cipal axes parallel to th e axes of the refere nce coordin a te sys te m , it will necessarily undergo a se ve re strain at th e boundaries of th e c ube in whic h th e particles are confined . It the refore see ms likely that homogeneous nucle ation has occurred in these run s but tha t th e failure to reach a perfect cr ystallin e structure is an artifact of the periodic boundary conditions and perhaps of the small size of the sys te m.
Similar computations on the T= 2.74 isoth erm we re performed at p = 1.30 which is approximately 20 percent beyond the freezing d ensity [2]. The res ults are give n in figures 3a and 3b. The gen eral be havior of the pressure and e nergy is simila r to the pre vious simulation. (The pressure and e ne rgy after conve rgence are 59 and -2.4, respectively, co mpared to 53 a nd -3.3 for a simulation of an fcc array at th e same te mperature and de nsity.) Th e molecular order, howe ve r , s uggests differe nt features. Th e pair correlation fun ction in fi gure 3b shows no e vide nce of the fcc orderin g indicated by th e solid c urve in fi gure 2b. Vis ual in s pection of the atomic coordin ates at 19 X 10 6 configurati ons by means of ste reoscopi c pairs of ph otographs indi cated an amorphou s stru cture, with little e vid e nce of di storted crys tallin e planes. During th e same sequence of confi gurations, th e molecular stru cture gradually ceased to be la bile a nd th e particles l curve of figure 2b. In previous studies [2] we have observed that the existence of fcc lattice structure results in a subdued peak or "step" in this region of the pair correlation function. The pair correlation shown in figure 3b exhibits some of the features found in studies of random close-packed arrays of hard spheres [3]. In particular, the broad second maximum can be interpreted as two separate peaks, centered at approximately r= 1.75 and r= 1.95. Bernal [4] attributes the first of these to the presence of tetrahedra which share a common base, and the second to the collineation of three particles. The relation between the structure of a random close-packed system and that of an ideal glass has been discussed in detail by Turnbull [5]. It is sufficient to remark here that the structure of the glass state in macroscopic systems is still an important unresolved problem.
We have investigated the reproducibility of the simulations in the metastable region by repeating the calculations with initial coordinates obtained trom a random number generator. The results for p = 1.

19
densities [2]. The dashed curves indicate the simulations are reproducible to within 4 percent and we conjecture that the differences would be smaller if the simulation in figure 4a had been continued to the same number of configurations as that in figure 3a. her:ame locked onto a disordered array. This was also manifested in the pair correlation function which did not change by more than ± 1 percent beyond 14 X 10 6 configurations.
We infer from this that the molecular order became locally more close-packed as the pressure and energy converged. The amorphous structure at convergence is suggestive of the formation of a glassy state. In figure 3b there is no evidence of the subdued second The pair correlation function obtained from the simulation in figure 4a at 15 X 10 6 is quite similar to that shown in figure 3b. In particular the function shows no evidence of the crystallinity shown in figures Ib and 2b.
The constant pressure heat capacity for the Lennard-Jones sys tem c an be computed in terms of press ure and energy fluctuation s, I1P and I1U respectively, by

«(I1P) 2) /kBT-7PfV+8U/VZ'
The first term on the ri ght is the con stant volum e heat capacity and the seco nd te rm re prese nts the ratio of the square of the compressibility to the expansivity. For potentials of the Lennard-Jones type (integer powe r law fun ction) th e ratio can be wri tte n in te rms of pressure and e nergy flu ctuations [8].
The results for th e co nfi gurational heat capacity for the simul a tion s hown in fi gure 4a are giv en in fi gure 4b. Since the press ure a nd energy flu ctuations are large, the heat capacity also has large flu ctuations. For reference, we note that th e heat capacity of the crystalline solid was es timated to be 3.4 whic h is close to the res ult fo r the glassy state s hown in fi gure 4b.