Communication: Phase diagram of C36 by atomistic molecular dynamics and thermodynamic integration through coexistence regions M. C. Abramo, C. Caccamo, D. Costa, and G. Munaò Citation: The Journal of Chemical Physics 141, 091103 (2014); doi: 10.1063/1.4894809 View online: http://dx.doi.org/10.1063/1.4894809 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/9?ver=pdfcov Published by the AIP Publishing Articles you may be interested in On the existence of a third-order phase transition beyond the Andrews critical point: A molecular dynamics study J. Chem. Phys. 135, 224506 (2011); 10.1063/1.3666848 Phase diagram of Hertzian spheres J. Chem. Phys. 131, 044514 (2009); 10.1063/1.3186742 Thermodynamic behavior and structural properties of an aqueous sodium chloride solution upon supercooling J. Chem. Phys. 128, 244508 (2008); 10.1063/1.2939118 Phase diagrams of alkali halides using two interaction models: A molecular dynamics and free energy study J. Chem. Phys. 126, 024503 (2007); 10.1063/1.2423030 Free energy determination of phase coexistence in model C 60 : A comprehensive Monte Carlo study J. Chem. Phys. 118, 304 (2003); 10.1063/1.1525800
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
THE JOURNAL OF CHEMICAL PHYSICS 141, 091103 (2014)
Communication: Phase diagram of C36 by atomistic molecular dynamics and thermodynamic integration through coexistence regions M. C. Abramo, C. Caccamo,a) D. Costa, and G. Munaò Dipartimento di Fisica e di Scienze della Terra, Università degli Studi di Messina and CNISM (Consorzio Nazionale Interuniversitario di Struttura della Materia) Viale F. Stagno d’Alcontres 31, 98166 Messina, Italy
(Received 10 July 2014; accepted 26 August 2014; published online 4 September 2014) We report an atomistic molecular dynamics determination of the phase diagram of a rigid-cage model of C36 . We first show that free energies obtained via thermodynamic integrations along isotherms displaying “van der Waals loops,” are fully reproduced by those obtained via isothermal-isochoric integration encompassing only stable states. We find that a similar result also holds for isochoric paths crossing van der Waals regions of the isotherms, and for integrations extending to rather high densities where liquid-solid coexistence can be expected to occur. On such a basis we are able to map the whole phase diagram of C36 , with resulting triple point and critical temperatures about 1770 K and 2370 K, respectively. We thus predict a 600 K window of existence of a stable liquid phase. Also, at the triple point density, we find that the structural functions and the diffusion coefficient maintain a liquid-like character down to 1400–1300 K, this indicating a wide region of possible supercooling. We discuss why all these features might render possible the observation of the melting of C36 fullerite and of its liquid state, at variance with what previously experienced for C60 . © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4894809] Atomistic simulations of rigid-cage models of fullerenes have been performed in the past, in particular for C60 and C70 crystalline fullerites, mostly in relation with solid state properties of these materials not far from ambient conditions.1–6 More recently two of us7 have reported the results of an atomistic Molecular Dynamics (MD) investigation of a rigid molecule model of C36 fullerite8, 9 in which both ambient and high-temperature properties have been addressed. Extensive simulations and theoretical studies of central potential models of C60 and C70 have also been reported (see Refs. 10–17 and references therein) but, as shown in Ref. 7, such a representation turns out unsuitable for C36 . In this Communication, we report the determination of the phase diagram of C36 , performed by means of atomistic MD in which the same model of Ref. 7 is employed: carbon atoms are fixed over a rigid molecular cage according to the experimental coordinates, to form the structure shown in Fig. 1; carbons on different cages are then assumed to interact via a Lennard-Jones potential whose parameters are fixed by imposing the fit7 of physical quantities of the solid C36 fullerite.8, 9 We calculate via thermodynamic integration (TI) the free energy necessary for the determination of phase coexistence lines. As is well known, numerical calculations of this type are rather demanding and, to the best of our knowledge, TI simulation calculations in the present context have been reported only for the hydration of fullerenes (see Refs. 18 and 19 and references therein). Actually, beside the computational effort, the TI procedure is usually considered hampered by the
a) Author to whom correspondence should be addressed. Electronic mail:
[email protected]
0021-9606/2014/141(9)/091103/5/$30.00
possible crossing of “metastable” or “unstable” states. However, recent simulation studies (see Refs. 20 and 21 and references therein) have enlightened the effective nature of such states in connection with “van der Waals loops” displayed by the chemical potentials, the occurrence of phase separation inside the simulation box, and the finite size of this latter. Here we report calculations based on isothermal paths exhibiting van der Waals-like loops, and on isothermal-isochoric (mixed) paths which run only over stable states. We show that the free energies calculated through the two routes are practically coincident, and that this also happens for mixed paths crossing the van der Waals region, or high-density regions where the liquid-solid coexistence usually occurs. We discuss the emerging equivalence of the free energies also in the light of the previously cited studies on metastable and unstable states in simulation of finite size systems.20, 21 Hinging on the above results, we determine both the liquid-vapor and the liquid-solid coexistence lines, and thereby get estimates of the critical and triple point parameters of C36 . It turns out that the triple and critical point temperatures, Ttr and Tcr , are significantly lower, and much higher, respectively, of their counterparts for a rigid cage model of C60 .10, 11, 13, 15–17 We get thus evidence of a huge expansion of the stable liquid phase pocket of C36 in comparison to that early found in C60 .10–12, 14, 16, 17 Moreover, the diffusion coefficient and the structural arrangement at the triple point density of C36 , maintain liquid-like features over several hundreds degrees below Ttr , by thus indicating a substantial extension of the “supercooled” liquid region. We recall in this instance that attempts to observe experimentally the liquid phase of C60 by heating its fullerite at high temperature,22 failed because of the amorphization of the C60 cage, occurring ∼700 K below the estimated Ttr .10, 16, 17 The reasons of such an instability of
141, 091103-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-2
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
1100 2500K
180
2350K
900
2200K −20
2000K
700
1800K
P[MPa]
−220
2
500
50
300
0
100
−50 0.0
2.2
0.4
2.4
0.8
1.2
2.6
1.6
1700K
2.8
2.0
2.4
−100
the C60 molecule in the solid phase have been object of debate in the past (see Ref. 23 and references therein). Here we find that liquid-like, albeit metastable, conditions may be encountered down to temperatures comparable to those for which the amorphization phenomenon has been observed; such a circumstance seems of interest as far as a possible detection of the C36 fullerite melting is concerned. We recall the Lennard-Jones potential parameters for the C–C interaction, = 90 K, in kB (Boltzmann constant) units, and σ = 0.2816 nm, employed in Ref. 7 (to be consulted for other details). NVT MD simulations for C36 have been performed with such parameters, by employing the MOLDY code.24 The simulation samples have been constituted with 1000 C36 rigid molecules enclosed in cubic boxes of edge cell length varying according to the number density ρ in the range 0.01 ≤ ρ ≤ 2.7 nm−3 (ρ = 3.086 nm−3 being the density of the solid at ambient conditions). An integration time step of 5 × 10−15 s has been adopted. All simulations have been started from the same initial configuration and a Nosé-Poincaré thermostat has been used.25 A preliminary NVE MD simulation has always been performed for 100 ps, followed by NVT MD runs for 1000 ps and, in some cases, up to 1500 ps, by cumulating averages over the last 250−500 ps (depending on density and temperature). Errors have been estimated through block averages each of 50 ps. The isotherms T = 2500, 2350, 2200, 2000, 1800, 1700, 1500 K have been investigated for 28 densities. The triple point isochore ρ = 2.2 nm−3 has also been studied for T = 1400, 1300, 1200 K. Pressures as a function of the density are shown in Fig. 2. As visible, all isotherms for T < 2500 K are constituted by two non-monotonic portions, having the appearance of van der Waals loops, in the density intervals 0 < ρ < 2.2 and 2.2 < ρ < 2.7 nm−3 . This is illustrated in detail in the insets of Fig. 2 for the T = 1800 K case. The T = 2500 K pattern also exhibits a loop at high densities while, in the low-density range (bottom inset), the curve is monotonic with a horizontal inflection point, as typical of a critical or marginally critical isotherm. The excess free energy per particle is then calculated ρ through an isothermal integration: βFex /N = 0 (βP /ρ −1 − 1)/ρ dρ where β = (kB T) . Results for various isotherms are shown in Fig. 3. Since the isotherms display a van-derWaals shape, one can plausibly assume that the related in-
−300 0.0
0.5
1.0
1.5 −3 ρ[nm ]
2.0
2.5
3.0
FIG. 2. Equation of state of C36 at various temperatures. (Insets) (top) isotherm at T = 1800 K and at high density on expanded scale; (bottom) isotherms at T = 2500 and 1800 K at low density on expanded scale. Error bars, displayed in the insets for two temperatures, have similar magnitude for the other T’s. Lines are guides to the eye.
tegration paths encompass metastable or unstable states, and that this might bias the reliability of the calculated βFex /N. We have then estimated this quantity also through a mixed isothermal-isochoric path, as follows: an integration is first performed along the T = 2500 K isotherm, for which no loops occur, up to the relatively high density ρ = 2.2 nm−3 ; then, βFex /N for T < 2500 K are determined for densities 1.8 ≤ ρ ≤ 2.2 nm−3 by adding to the 2500 K value the isochoric inteβ gration contribution Fex = β U (β ) dβ where U is the ex0
cess internal energy. Here β 0 = (kB T0 )−1 and T0 = 2500 K. As it will appear below, the outlined procedure, applied over the above density-temperature range, actually allows us to reach thermodynamic points located just beyond the binodal line without crossing this latter. As visible in Fig. 3, the βFex /N
0
−1
βFex/N
FIG. 1. The C36 model employed in this work. The molecular size is ∼0.5 nm.
−2
−3
−4 0.0
0.5
1.0
1.5
2.0
2.5
3.0
−3
ρ[nm ] FIG. 3. Comparison of excess free energies as obtained through isothermal (circles) and isothermal-isochoric (squares) integrations, at T = 2350, 2200, 2000, 1800 K, from top to bottom.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
0
−1
2000K 1700K
βμ
−2
−3
−4
−5 −100
0
100
200
P[MPa] FIG. 4. μ vs P at constant T. Diamonds and triangles mark the intersections of the μ branches corresponding to liquid-vapor and liquid-solid coexistence, respectively. Errors on μ are of the same entity as those on βFex /N.
2600
2200 500
1800
P[MPa]
calculated through the two different paths practically overlap to within the error. We then extend the isochoric integration to intermediate densities 0.5 ≤ ρ ≤ 1.8 nm−3 by going significantly below 2500 K. As shown, also in this case we find that the isothermally determined βFex /N are well reproduced. Finally, the T = 2500 K integration is extended to the 2.2 ≤ ρ ≤ 2.7 nm−3 range, thus crossing the non-monotonic portion of this isotherm at high densities; isochoric integrations are then performed at various densities down to the different temperatures explored. The βFex /N thereby obtained again compare quite well with their isothermal counterparts, as visible in Fig. 3. We interpret the reproducibility of βFex /N calculated through the different paths, as an evidence that integrations for T < 2500 K and up to high densities actually extend over stable thermodynamic states. This conclusion looks consistent with what pointed out in simulations of Refs. 20 and 21, where the authors document that although van der Waals-like loops are visible in the chemical potential, the system actually samples thermodynamically stable states characterized by phase separation occurring inside the box because of finite size effects. Hinging on the estimates of βFex /N we construct the chemical potential βμ = βF/N + βP/ρ, where the Helmholtz free energy βF/N is obtained from βFex /N by adding the ideal contribution. The βμ vs P patterns are shown in Fig. 4, limited to two temperatures for clarity purposes. As visible, two distinct crossings in the μ pattern are present at each temperature, yielding two equilibrium pressures Peq to be associated with liquid-vapor (LV) and liquid-solid (LS) coexistence. From the Peq the coexisting densities can be established according to Fig. 2; the resulting binodal and freezing-melting loci are displayed in Fig. 5. Note that because of the flat intersection in the μ vs P pattern at the LS equilibrium, the error in Peq is larger than for the LV case and this reflects in the coexisting densities error. The locus of LV coexistence points appears slightly asymmetric toward the low density end, while the melting
T[K]
091103-3
300
100
−100
1600
2100
2600
T[K]
1400 0.0
0.5
1.0
1.5 −3 ρ[nm ]
2.0
2.5
3.0
FIG. 5. Phase diagram of C36 in the T − ρ representation. Full line: fit of the MD binodal points (circles) described in the text. Freezing: squares; melting: diamonds; critical and triple points: up and down triangles, respectively. Error bars, if not displayed, are smaller than the size of the symbols. (Inset) P vs T phase coexistence: liquid-vapor, circles; liquid-solid, diamonds.
and freezing loci have fairly usual appearance. We estimate Ttr = 1770 ± 20 K with ρ tr = 2.26 ± 0.02 nm−3 , and Ptr = 6 ± 3 MPa. As for the critical parameters, we perform an interpolation of the available coexisting liquid-vapor densities, ρ L and ρV , by means of the scaling law ρL − ρV ∝ |T − Tcr |β , wherein the non-classical critical exponent β = 0.32 is employed, and of the law of rectilinear diameters.28 The resulting fit is displayed in the figure. We obtain Tcr 2370 K and ρ cr 0.92 nm−3 . We estimate on this basis a critical pressure Pcr = 10 ± 0.4 MPa. It follows that the range of existence of a thermodynamically stable liquid phase has an amplitude of ∼600 K, much more indeed than the 80 K interval found for the C60 case.10, 13, 15–17 Since the C–C interaction in the present case has an attractive well depth = 90 K, almost three times greater than the one adopted in C60 (see Refs. 7 and 16), a certain increase of Tcr in C36 might be expected; the magnitude of the overall expansion of the liquid pocket appears nonetheless remarkable. The structural characterization of the different phases is shown in Fig. 6. The center-center radial distribution functions g(r) are practically zero for r < 0.5 nm, consistently with the size of the rigid molecular cage (see Fig. 1). They also display fully liquid-like features for thermodynamic points inside the liquid pocket (e.g., T = 1800 K, ρ = 2.2 nm−3 ). Close to freezing at constant T (panel (a)), a shoulder begins to develop in the second maximum of g(r), as visible for ρ 2.45 nm−3 at r = 1.35 nm, and well defined peaks, typical of crystalline ordering, emerge in the approach to, and crossing of, the melting locus at ρ 2.6 nm−3 . Also, if one follows the structural evolution for T < Ttr at fixed ρ = ρ tr (panel (b)) it appears that liquid-like configurations persist down to 1400 K, this signaling a wide extension of the liquid supercooling region. We shall further comment on this point when discussing dynamical features. The analysis of microscopic configurations at T = 1800 K and different densities (not shown here for space reasons) documents the presence of a vapor at low densities,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-4
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
6 (a)
2.200 2.449 2.519 2.591
(b)
1500 1400 1300 1200
4
g(r)
2 0 10 8 6 4 2 0
0
1
2
3
4
r[nm]
FIG. 6. Center-center radial distribution function g(r). Panel (a) constant temperature T = 1800 K and different densities (in nm−3 ). Panel (b) constant density ρ = 2.2 nm−3 and different temperatures (in K).
of a liquid at ρ = 2.2 nm−3 , and the onset of a crystalline arrangement at ρ = 2.591 nm−3 , respectively. At intermediate densities, ρ = 1.5 nm−3 say, the system displays a nonhomogeneous arrangement with the coexistence of a liquid and a vapor region without a well defined interface. The dynamical characterization of the liquid and solid phases can then be appreciated in Fig. 7, where the mean square displacement vs the simulation time is displayed. We estimate from such plots the diffusion coefficient D: inside the liquid pocket and close to the triple point we find D 2 × 10−5 cm2 /s, a value typical of simple liquids. Upon increasing the density and approaching the melting line, as visible in panel (a) of Fig. 7, one finds that D undergoes a drop of two orders of magnitude as compatible with the onset of the solid arrangement. One also finds, for the sequence visible in panel (b) of Fig. 7, that D 10−5 cm2 /s down to T 1500 K, and
that it attains a fraction of this value at 1400 K. In agreement with the structural evidences before discussed, this signals a persistence of supercooling at least ∼350 K below Ttr . This circumstance might reveal precious as for the possible observation of C36 in the liquid state, should any amorphization of the C36 heated fullerite occur at the same temperature as in C60 , i.e., at ∼1300 K, competitive with, or possibly higher than the lower temperature boundary of liquid-like behavior. We observe that in relation to the LS coexistence line, an assessment of the βFex /N estimate inside the solid phase might be performed in terms, for instance, of TI starting from an Einstein solid estimate of the free energy, or of some theoretical approach, as done in previous works for C60 .16, 17 Therein, however, the fullerene was described in terms of a central potential model, and the extension of those methods to our atomistic C36 does not appear of immediate feasibility. A similar remark can be made also for the determination of the freezing line in terms of the multiparticle residual entropy behavior,26 an approach exploited by us in previous works for a central potential model of C60 and for a simple fluid, respectively.16, 27 The binodal line could also be derived via direct methods as for instance successive umbrella sampling, grand canonical, or Gibbs ensemble Monte Carlo simulations.28 On the other hand, the advantage of the present MD study is to provide a full characterization of the fluid phase, including the freezing and melting lines, as well as the dynamical properties of the system. We plan to apply the calculation of the phase diagram as here described to other multi-site molecules as, for instance, mesoscopic models of protein solutions recently investigated by us.29, 30 We warmly thank Professor Paolo V. Giaquinta for stimulating discussions and a critical reading of the paper. Useful contacts with Professor Giovanni Ciccotti are acknowledged. Support from the PRIN-MIUR 2010–2011 project is also acknowledged. C.C. acknowledges Grant No. 91796 (Blue Skies) from the National Research Foundation (NRF) of South Africa.
2 2.200 2.318 2.383 2.449 2.519 2.591
1 A.
0 6 1800 1500 1400 1300 1200
4
(b)
2 0
0
50
Cheng and M. L. Klein, J. Phys. Chem. 95, 6750 (1991). Cheng and M. L. Klein, Phys. Rev. B 45, 1889 (1992). 3 M. Sprik, A. Cheng, and M. L. Klein, J. Phys. Chem. 96, 2027 (1992). 4 M. Sprik, A. Cheng, and M. L. Klein, Phys. Rev. Lett. 69, 1660 (1992). 5 M. Sprik and M. L. Klein, J. Phys. Chem. 98, 9297 (1994). 6 M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane, and R. Ruberto, Phys. Rev. E 69, 031112 (2004). 7 M. C. Abramo and C. Caccamo, J. Chem. Phys. 128, 074503 (2008). 8 C. Piskoti, T. Yarger, and A. Zettl, Nature (London) 393, 771 (1998). 9 M. Côté, J. C. Grossman, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 81, 697 (1998). 10 A. Cheng, M. L. Klein, and C. Caccamo, Phys. Rev. Lett. 71, 1200 (1993). 11 M. H. J. Hagen, E. J. Meijer, G. C. A. M. Mooij, D. Frenkel, and H. N. W. Lekkerkerker, Nature (London) 365, 425 (1993). 12 N. W. Ashcroft, Nature (London) 365, 387 (1993). 13 C. Caccamo, Phys. Rev. B 51, 3387 (1995). 14 M. Hasegawa and K. Ohno, J. Chem. Phys. 111, 5955 (1999). 15 M. C. Abramo, C. Caccamo, D. Costa, and G. Pellicane, Europhys. Lett. 54, 468 (2001). 16 D. Costa, G. Pellicane, C. Caccamo, and M. C. Abramo, J. Chem. Phys. 118, 304 (2003). 2 A.
2
MSD[nm ]
1
(a)
100 150 time[ps]
200
250
FIG. 7. Mean square displacement vs simulation time. Panel (a) constant T = 1800 K and different densities (in nm−3 ). Panel (b) constant ρ = 2.2 nm−3 and different temperatures (in K).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-5 17 D.
Abramo et al.
Costa, G. Pellicane, C. Caccamo, E. Schöll-Paschinger, and G. Kahl, Phys. Rev. E 68, 021104 (2003). 18 C. Maciel, E. E. Fileti, and R. Rivelino, Chem. Phys. Lett. 507, 244 (2011). 19 G. Colherinhas, T. L. Fonseca, and E. E. Fileti, Carbon 49, 187 (2011). 20 K. Binder, B. J. Block, P. Virnau, and A. Tröster, Am. J. Phys. 80, 1099 (2012). 21 B. G. MacDowell, P. Virnau, M. Müller, and K. Binder, J. Chem. Phys. 120, 5293 (2004). 22 M. R. Stetzer, P. A. Heiney, J. E. Fischer, and A. R. McGhie, Phys. Rev. B 55, 127 (1997). 23 M. C. Abramo and C. Caccamo, J. Chem. Phys. 106, 6475 (1997).
J. Chem. Phys. 141, 091103 (2014) 24 K.
Refson, Comput. Phys. Commun. 126, 310 (2000). D. Bond, B. J. Leimkuhler, and B. B. Laird, J. Comput. Phys. 151, 114 (1999). 26 P. V. Giaquinta and G. Giunta, Physica A 187, 145 (1992). 27 C. Caccamo, D. Costa, and G. Pellicane, J. Chem. Phys. 109, 4498 (1998). 28 D. Frenkel and B. Smith, Understanding Molecular Simulation (Academic, London, 1996). 29 F. Carlsson, P. Linse, and M. Malmsten, J. Phys. Chem. B 105, 9040 (2001). 30 M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane, and R. Ruberto, J. Phys. Chem. B 114, 9109 (2010). 25 S.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
THE JOURNAL OF CHEMICAL PHYSICS 141, 091103 (2014)
Communication: Phase diagram of C36 by atomistic molecular dynamics and thermodynamic integration through coexistence regions M. C. Abramo, C. Caccamo,a) D. Costa, and G. Munaò Dipartimento di Fisica e di Scienze della Terra, Università degli Studi di Messina and CNISM (Consorzio Nazionale Interuniversitario di Struttura della Materia) Viale F. Stagno d’Alcontres 31, 98166 Messina, Italy
(Received 10 July 2014; accepted 26 August 2014; published online 4 September 2014) We report an atomistic molecular dynamics determination of the phase diagram of a rigid-cage model of C36 . We first show that free energies obtained via thermodynamic integrations along isotherms displaying “van der Waals loops,” are fully reproduced by those obtained via isothermal-isochoric integration encompassing only stable states. We find that a similar result also holds for isochoric paths crossing van der Waals regions of the isotherms, and for integrations extending to rather high densities where liquid-solid coexistence can be expected to occur. On such a basis we are able to map the whole phase diagram of C36 , with resulting triple point and critical temperatures about 1770 K and 2370 K, respectively. We thus predict a 600 K window of existence of a stable liquid phase. Also, at the triple point density, we find that the structural functions and the diffusion coefficient maintain a liquid-like character down to 1400–1300 K, this indicating a wide region of possible supercooling. We discuss why all these features might render possible the observation of the melting of C36 fullerite and of its liquid state, at variance with what previously experienced for C60 . © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4894809] Atomistic simulations of rigid-cage models of fullerenes have been performed in the past, in particular for C60 and C70 crystalline fullerites, mostly in relation with solid state properties of these materials not far from ambient conditions.1–6 More recently two of us7 have reported the results of an atomistic Molecular Dynamics (MD) investigation of a rigid molecule model of C36 fullerite8, 9 in which both ambient and high-temperature properties have been addressed. Extensive simulations and theoretical studies of central potential models of C60 and C70 have also been reported (see Refs. 10–17 and references therein) but, as shown in Ref. 7, such a representation turns out unsuitable for C36 . In this Communication, we report the determination of the phase diagram of C36 , performed by means of atomistic MD in which the same model of Ref. 7 is employed: carbon atoms are fixed over a rigid molecular cage according to the experimental coordinates, to form the structure shown in Fig. 1; carbons on different cages are then assumed to interact via a Lennard-Jones potential whose parameters are fixed by imposing the fit7 of physical quantities of the solid C36 fullerite.8, 9 We calculate via thermodynamic integration (TI) the free energy necessary for the determination of phase coexistence lines. As is well known, numerical calculations of this type are rather demanding and, to the best of our knowledge, TI simulation calculations in the present context have been reported only for the hydration of fullerenes (see Refs. 18 and 19 and references therein). Actually, beside the computational effort, the TI procedure is usually considered hampered by the
a) Author to whom correspondence should be addressed. Electronic mail:
[email protected]
0021-9606/2014/141(9)/091103/5/$30.00
possible crossing of “metastable” or “unstable” states. However, recent simulation studies (see Refs. 20 and 21 and references therein) have enlightened the effective nature of such states in connection with “van der Waals loops” displayed by the chemical potentials, the occurrence of phase separation inside the simulation box, and the finite size of this latter. Here we report calculations based on isothermal paths exhibiting van der Waals-like loops, and on isothermal-isochoric (mixed) paths which run only over stable states. We show that the free energies calculated through the two routes are practically coincident, and that this also happens for mixed paths crossing the van der Waals region, or high-density regions where the liquid-solid coexistence usually occurs. We discuss the emerging equivalence of the free energies also in the light of the previously cited studies on metastable and unstable states in simulation of finite size systems.20, 21 Hinging on the above results, we determine both the liquid-vapor and the liquid-solid coexistence lines, and thereby get estimates of the critical and triple point parameters of C36 . It turns out that the triple and critical point temperatures, Ttr and Tcr , are significantly lower, and much higher, respectively, of their counterparts for a rigid cage model of C60 .10, 11, 13, 15–17 We get thus evidence of a huge expansion of the stable liquid phase pocket of C36 in comparison to that early found in C60 .10–12, 14, 16, 17 Moreover, the diffusion coefficient and the structural arrangement at the triple point density of C36 , maintain liquid-like features over several hundreds degrees below Ttr , by thus indicating a substantial extension of the “supercooled” liquid region. We recall in this instance that attempts to observe experimentally the liquid phase of C60 by heating its fullerite at high temperature,22 failed because of the amorphization of the C60 cage, occurring ∼700 K below the estimated Ttr .10, 16, 17 The reasons of such an instability of
141, 091103-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-2
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
1100 2500K
180
2350K
900
2200K −20
2000K
700
1800K
P[MPa]
−220
2
500
50
300
0
100
−50 0.0
2.2
0.4
2.4
0.8
1.2
2.6
1.6
1700K
2.8
2.0
2.4
−100
the C60 molecule in the solid phase have been object of debate in the past (see Ref. 23 and references therein). Here we find that liquid-like, albeit metastable, conditions may be encountered down to temperatures comparable to those for which the amorphization phenomenon has been observed; such a circumstance seems of interest as far as a possible detection of the C36 fullerite melting is concerned. We recall the Lennard-Jones potential parameters for the C–C interaction, = 90 K, in kB (Boltzmann constant) units, and σ = 0.2816 nm, employed in Ref. 7 (to be consulted for other details). NVT MD simulations for C36 have been performed with such parameters, by employing the MOLDY code.24 The simulation samples have been constituted with 1000 C36 rigid molecules enclosed in cubic boxes of edge cell length varying according to the number density ρ in the range 0.01 ≤ ρ ≤ 2.7 nm−3 (ρ = 3.086 nm−3 being the density of the solid at ambient conditions). An integration time step of 5 × 10−15 s has been adopted. All simulations have been started from the same initial configuration and a Nosé-Poincaré thermostat has been used.25 A preliminary NVE MD simulation has always been performed for 100 ps, followed by NVT MD runs for 1000 ps and, in some cases, up to 1500 ps, by cumulating averages over the last 250−500 ps (depending on density and temperature). Errors have been estimated through block averages each of 50 ps. The isotherms T = 2500, 2350, 2200, 2000, 1800, 1700, 1500 K have been investigated for 28 densities. The triple point isochore ρ = 2.2 nm−3 has also been studied for T = 1400, 1300, 1200 K. Pressures as a function of the density are shown in Fig. 2. As visible, all isotherms for T < 2500 K are constituted by two non-monotonic portions, having the appearance of van der Waals loops, in the density intervals 0 < ρ < 2.2 and 2.2 < ρ < 2.7 nm−3 . This is illustrated in detail in the insets of Fig. 2 for the T = 1800 K case. The T = 2500 K pattern also exhibits a loop at high densities while, in the low-density range (bottom inset), the curve is monotonic with a horizontal inflection point, as typical of a critical or marginally critical isotherm. The excess free energy per particle is then calculated ρ through an isothermal integration: βFex /N = 0 (βP /ρ −1 − 1)/ρ dρ where β = (kB T) . Results for various isotherms are shown in Fig. 3. Since the isotherms display a van-derWaals shape, one can plausibly assume that the related in-
−300 0.0
0.5
1.0
1.5 −3 ρ[nm ]
2.0
2.5
3.0
FIG. 2. Equation of state of C36 at various temperatures. (Insets) (top) isotherm at T = 1800 K and at high density on expanded scale; (bottom) isotherms at T = 2500 and 1800 K at low density on expanded scale. Error bars, displayed in the insets for two temperatures, have similar magnitude for the other T’s. Lines are guides to the eye.
tegration paths encompass metastable or unstable states, and that this might bias the reliability of the calculated βFex /N. We have then estimated this quantity also through a mixed isothermal-isochoric path, as follows: an integration is first performed along the T = 2500 K isotherm, for which no loops occur, up to the relatively high density ρ = 2.2 nm−3 ; then, βFex /N for T < 2500 K are determined for densities 1.8 ≤ ρ ≤ 2.2 nm−3 by adding to the 2500 K value the isochoric inteβ gration contribution Fex = β U (β ) dβ where U is the ex0
cess internal energy. Here β 0 = (kB T0 )−1 and T0 = 2500 K. As it will appear below, the outlined procedure, applied over the above density-temperature range, actually allows us to reach thermodynamic points located just beyond the binodal line without crossing this latter. As visible in Fig. 3, the βFex /N
0
−1
βFex/N
FIG. 1. The C36 model employed in this work. The molecular size is ∼0.5 nm.
−2
−3
−4 0.0
0.5
1.0
1.5
2.0
2.5
3.0
−3
ρ[nm ] FIG. 3. Comparison of excess free energies as obtained through isothermal (circles) and isothermal-isochoric (squares) integrations, at T = 2350, 2200, 2000, 1800 K, from top to bottom.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
0
−1
2000K 1700K
βμ
−2
−3
−4
−5 −100
0
100
200
P[MPa] FIG. 4. μ vs P at constant T. Diamonds and triangles mark the intersections of the μ branches corresponding to liquid-vapor and liquid-solid coexistence, respectively. Errors on μ are of the same entity as those on βFex /N.
2600
2200 500
1800
P[MPa]
calculated through the two different paths practically overlap to within the error. We then extend the isochoric integration to intermediate densities 0.5 ≤ ρ ≤ 1.8 nm−3 by going significantly below 2500 K. As shown, also in this case we find that the isothermally determined βFex /N are well reproduced. Finally, the T = 2500 K integration is extended to the 2.2 ≤ ρ ≤ 2.7 nm−3 range, thus crossing the non-monotonic portion of this isotherm at high densities; isochoric integrations are then performed at various densities down to the different temperatures explored. The βFex /N thereby obtained again compare quite well with their isothermal counterparts, as visible in Fig. 3. We interpret the reproducibility of βFex /N calculated through the different paths, as an evidence that integrations for T < 2500 K and up to high densities actually extend over stable thermodynamic states. This conclusion looks consistent with what pointed out in simulations of Refs. 20 and 21, where the authors document that although van der Waals-like loops are visible in the chemical potential, the system actually samples thermodynamically stable states characterized by phase separation occurring inside the box because of finite size effects. Hinging on the estimates of βFex /N we construct the chemical potential βμ = βF/N + βP/ρ, where the Helmholtz free energy βF/N is obtained from βFex /N by adding the ideal contribution. The βμ vs P patterns are shown in Fig. 4, limited to two temperatures for clarity purposes. As visible, two distinct crossings in the μ pattern are present at each temperature, yielding two equilibrium pressures Peq to be associated with liquid-vapor (LV) and liquid-solid (LS) coexistence. From the Peq the coexisting densities can be established according to Fig. 2; the resulting binodal and freezing-melting loci are displayed in Fig. 5. Note that because of the flat intersection in the μ vs P pattern at the LS equilibrium, the error in Peq is larger than for the LV case and this reflects in the coexisting densities error. The locus of LV coexistence points appears slightly asymmetric toward the low density end, while the melting
T[K]
091103-3
300
100
−100
1600
2100
2600
T[K]
1400 0.0
0.5
1.0
1.5 −3 ρ[nm ]
2.0
2.5
3.0
FIG. 5. Phase diagram of C36 in the T − ρ representation. Full line: fit of the MD binodal points (circles) described in the text. Freezing: squares; melting: diamonds; critical and triple points: up and down triangles, respectively. Error bars, if not displayed, are smaller than the size of the symbols. (Inset) P vs T phase coexistence: liquid-vapor, circles; liquid-solid, diamonds.
and freezing loci have fairly usual appearance. We estimate Ttr = 1770 ± 20 K with ρ tr = 2.26 ± 0.02 nm−3 , and Ptr = 6 ± 3 MPa. As for the critical parameters, we perform an interpolation of the available coexisting liquid-vapor densities, ρ L and ρV , by means of the scaling law ρL − ρV ∝ |T − Tcr |β , wherein the non-classical critical exponent β = 0.32 is employed, and of the law of rectilinear diameters.28 The resulting fit is displayed in the figure. We obtain Tcr 2370 K and ρ cr 0.92 nm−3 . We estimate on this basis a critical pressure Pcr = 10 ± 0.4 MPa. It follows that the range of existence of a thermodynamically stable liquid phase has an amplitude of ∼600 K, much more indeed than the 80 K interval found for the C60 case.10, 13, 15–17 Since the C–C interaction in the present case has an attractive well depth = 90 K, almost three times greater than the one adopted in C60 (see Refs. 7 and 16), a certain increase of Tcr in C36 might be expected; the magnitude of the overall expansion of the liquid pocket appears nonetheless remarkable. The structural characterization of the different phases is shown in Fig. 6. The center-center radial distribution functions g(r) are practically zero for r < 0.5 nm, consistently with the size of the rigid molecular cage (see Fig. 1). They also display fully liquid-like features for thermodynamic points inside the liquid pocket (e.g., T = 1800 K, ρ = 2.2 nm−3 ). Close to freezing at constant T (panel (a)), a shoulder begins to develop in the second maximum of g(r), as visible for ρ 2.45 nm−3 at r = 1.35 nm, and well defined peaks, typical of crystalline ordering, emerge in the approach to, and crossing of, the melting locus at ρ 2.6 nm−3 . Also, if one follows the structural evolution for T < Ttr at fixed ρ = ρ tr (panel (b)) it appears that liquid-like configurations persist down to 1400 K, this signaling a wide extension of the liquid supercooling region. We shall further comment on this point when discussing dynamical features. The analysis of microscopic configurations at T = 1800 K and different densities (not shown here for space reasons) documents the presence of a vapor at low densities,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-4
Abramo et al.
J. Chem. Phys. 141, 091103 (2014)
6 (a)
2.200 2.449 2.519 2.591
(b)
1500 1400 1300 1200
4
g(r)
2 0 10 8 6 4 2 0
0
1
2
3
4
r[nm]
FIG. 6. Center-center radial distribution function g(r). Panel (a) constant temperature T = 1800 K and different densities (in nm−3 ). Panel (b) constant density ρ = 2.2 nm−3 and different temperatures (in K).
of a liquid at ρ = 2.2 nm−3 , and the onset of a crystalline arrangement at ρ = 2.591 nm−3 , respectively. At intermediate densities, ρ = 1.5 nm−3 say, the system displays a nonhomogeneous arrangement with the coexistence of a liquid and a vapor region without a well defined interface. The dynamical characterization of the liquid and solid phases can then be appreciated in Fig. 7, where the mean square displacement vs the simulation time is displayed. We estimate from such plots the diffusion coefficient D: inside the liquid pocket and close to the triple point we find D 2 × 10−5 cm2 /s, a value typical of simple liquids. Upon increasing the density and approaching the melting line, as visible in panel (a) of Fig. 7, one finds that D undergoes a drop of two orders of magnitude as compatible with the onset of the solid arrangement. One also finds, for the sequence visible in panel (b) of Fig. 7, that D 10−5 cm2 /s down to T 1500 K, and
that it attains a fraction of this value at 1400 K. In agreement with the structural evidences before discussed, this signals a persistence of supercooling at least ∼350 K below Ttr . This circumstance might reveal precious as for the possible observation of C36 in the liquid state, should any amorphization of the C36 heated fullerite occur at the same temperature as in C60 , i.e., at ∼1300 K, competitive with, or possibly higher than the lower temperature boundary of liquid-like behavior. We observe that in relation to the LS coexistence line, an assessment of the βFex /N estimate inside the solid phase might be performed in terms, for instance, of TI starting from an Einstein solid estimate of the free energy, or of some theoretical approach, as done in previous works for C60 .16, 17 Therein, however, the fullerene was described in terms of a central potential model, and the extension of those methods to our atomistic C36 does not appear of immediate feasibility. A similar remark can be made also for the determination of the freezing line in terms of the multiparticle residual entropy behavior,26 an approach exploited by us in previous works for a central potential model of C60 and for a simple fluid, respectively.16, 27 The binodal line could also be derived via direct methods as for instance successive umbrella sampling, grand canonical, or Gibbs ensemble Monte Carlo simulations.28 On the other hand, the advantage of the present MD study is to provide a full characterization of the fluid phase, including the freezing and melting lines, as well as the dynamical properties of the system. We plan to apply the calculation of the phase diagram as here described to other multi-site molecules as, for instance, mesoscopic models of protein solutions recently investigated by us.29, 30 We warmly thank Professor Paolo V. Giaquinta for stimulating discussions and a critical reading of the paper. Useful contacts with Professor Giovanni Ciccotti are acknowledged. Support from the PRIN-MIUR 2010–2011 project is also acknowledged. C.C. acknowledges Grant No. 91796 (Blue Skies) from the National Research Foundation (NRF) of South Africa.
2 2.200 2.318 2.383 2.449 2.519 2.591
1 A.
0 6 1800 1500 1400 1300 1200
4
(b)
2 0
0
50
Cheng and M. L. Klein, J. Phys. Chem. 95, 6750 (1991). Cheng and M. L. Klein, Phys. Rev. B 45, 1889 (1992). 3 M. Sprik, A. Cheng, and M. L. Klein, J. Phys. Chem. 96, 2027 (1992). 4 M. Sprik, A. Cheng, and M. L. Klein, Phys. Rev. Lett. 69, 1660 (1992). 5 M. Sprik and M. L. Klein, J. Phys. Chem. 98, 9297 (1994). 6 M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane, and R. Ruberto, Phys. Rev. E 69, 031112 (2004). 7 M. C. Abramo and C. Caccamo, J. Chem. Phys. 128, 074503 (2008). 8 C. Piskoti, T. Yarger, and A. Zettl, Nature (London) 393, 771 (1998). 9 M. Côté, J. C. Grossman, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 81, 697 (1998). 10 A. Cheng, M. L. Klein, and C. Caccamo, Phys. Rev. Lett. 71, 1200 (1993). 11 M. H. J. Hagen, E. J. Meijer, G. C. A. M. Mooij, D. Frenkel, and H. N. W. Lekkerkerker, Nature (London) 365, 425 (1993). 12 N. W. Ashcroft, Nature (London) 365, 387 (1993). 13 C. Caccamo, Phys. Rev. B 51, 3387 (1995). 14 M. Hasegawa and K. Ohno, J. Chem. Phys. 111, 5955 (1999). 15 M. C. Abramo, C. Caccamo, D. Costa, and G. Pellicane, Europhys. Lett. 54, 468 (2001). 16 D. Costa, G. Pellicane, C. Caccamo, and M. C. Abramo, J. Chem. Phys. 118, 304 (2003). 2 A.
2
MSD[nm ]
1
(a)
100 150 time[ps]
200
250
FIG. 7. Mean square displacement vs simulation time. Panel (a) constant T = 1800 K and different densities (in nm−3 ). Panel (b) constant ρ = 2.2 nm−3 and different temperatures (in K).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
091103-5 17 D.
Abramo et al.
Costa, G. Pellicane, C. Caccamo, E. Schöll-Paschinger, and G. Kahl, Phys. Rev. E 68, 021104 (2003). 18 C. Maciel, E. E. Fileti, and R. Rivelino, Chem. Phys. Lett. 507, 244 (2011). 19 G. Colherinhas, T. L. Fonseca, and E. E. Fileti, Carbon 49, 187 (2011). 20 K. Binder, B. J. Block, P. Virnau, and A. Tröster, Am. J. Phys. 80, 1099 (2012). 21 B. G. MacDowell, P. Virnau, M. Müller, and K. Binder, J. Chem. Phys. 120, 5293 (2004). 22 M. R. Stetzer, P. A. Heiney, J. E. Fischer, and A. R. McGhie, Phys. Rev. B 55, 127 (1997). 23 M. C. Abramo and C. Caccamo, J. Chem. Phys. 106, 6475 (1997).
J. Chem. Phys. 141, 091103 (2014) 24 K.
Refson, Comput. Phys. Commun. 126, 310 (2000). D. Bond, B. J. Leimkuhler, and B. B. Laird, J. Comput. Phys. 151, 114 (1999). 26 P. V. Giaquinta and G. Giunta, Physica A 187, 145 (1992). 27 C. Caccamo, D. Costa, and G. Pellicane, J. Chem. Phys. 109, 4498 (1998). 28 D. Frenkel and B. Smith, Understanding Molecular Simulation (Academic, London, 1996). 29 F. Carlsson, P. Linse, and M. Malmsten, J. Phys. Chem. B 105, 9040 (2001). 30 M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane, and R. Ruberto, J. Phys. Chem. B 114, 9109 (2010). 25 S.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 12:55:23
