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.

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

FIG. 1.

The C36 model employed in this work. The molecular size is ∼0.5 nm.

FIG. 1.

The C36 model employed in this work. The molecular size is ∼0.5 nm.

Close modal

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,19 and references therein). Actually, beside the computational effort, the TI procedure is usually considered hampered by the possible crossing of “metastable” or “unstable” states. However, recent simulation studies (see Refs. 20,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 belowTtr, 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 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.

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.

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.

Close modal

The excess free energy per particle is then calculated through an isothermal integration:

$\beta F_{\rm ex}/N= \int _0^{\rho } (\beta P/\rho ^{\prime }\break -1)/\rho ^{\prime } \mathop {\rm \strut d\!}\nolimits \rho ^{\prime }$
βF ex /N=0ρ(βP/ρ1)/ρdρ where β = (kBT)−1. Results for various isotherms are shown in Fig. 3. Since the isotherms display a van-der-Waals shape, one can plausibly assume that the related integration 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 integration contribution
$\Delta F_{\rm ex}= \int _{\beta _0}^{\beta } U(\beta ^{\prime }) \mathop {\rm \strut d\!}\nolimits \beta ^{\prime }$
ΔF ex =β0βU(β)dβ
where U is the excess internal energy. Here β0 = (kBT0)−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 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.

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.

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.

Close modal

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,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.

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.

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.

Close modal
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.

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.

Close modal

The locus of LV coexistence points appears slightly asymmetric toward the low density end, while the melting 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∝|TTcr|β, 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,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.

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).

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).

Close modal

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, 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 non-homogeneous 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.

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).

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).

Close modal

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.

1.
A.
Cheng
and
M. L.
Klein
,
J. Phys. Chem.
95
,
6750
(
1991
).
2.
A.
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
).
17.
D.
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
).
24.
K.
Refson
,
Comput. Phys. Commun.
126
,
310
(
2000
).
25.
S. 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
).