CO₂ packing polymorphism under confinement in cylindrical nanopores

Abstract

We investigate the effect of cylindrical nano-confinement on the phase behaviour of a rigid model of carbon dioxide using both molecular dynamics and well tempered metadynamics. To this aim we study a simplified pore model across a parameter space comprising pore diameter, CO₂-pore wall potential and CO₂ density. In order to systematically identify ordering events within the pore model we devise a generally applicable approach based on the analysis of the distribution of intermolecular orientations. Our simulations suggest that, while confinement in nano-pores inhibits the formation of known crystal structures, it induces a remarkable variety of ordered packings unrelated to their bulk counterparts, and favours the establishment of short range order in the fluid phase. We summarise our findings by proposing a qualitative phase diagram for this model.

---

Design, System, Application

Understanding the effect of confinement on polymorph selection is key to direct the assembly of the next generation of self assembled molecular materials. In particular predicting the effect of confinement in inhibiting or enhancing the nucleation of ordered phases is pivotal to develop efficient crystallization processes. In this paper, using a combination of molecular dynamics and well tempered metadynamics, we analyse the effect of cylindrical nano confinement on the phase behaviour of CO₂. Our analysis is relevant to gain insight into the behaviour of CO₂ porous matrices with potential impact on carbon capture processes. Furthermore it allows to identify general principles of confinement-induced polymorph selection. For instance we find that confinement inhibit the nucleation of bulk crystal forms, while promoting the formation of distinct ordered packings. This work paves the way for a systematic exploration of confinement as a tool for directing the assembly of molecular materials.

1 Introduction

Confinement is known to play a role in the phase behaviour of molecular solids, most notably affecting polymorph selection.^1,2 For instance, a paradigmatic example of the dramatic effects of confinement on the spatial arrangement of molecules is provided by water, which, as proven both experimentally and computationally, displays a counter-intuitively complex phase diagram under confinement.^3–10 Understanding polymorphism in confined volumes is relevant both to describe natural processes^11,12 as well as for driving rational materials and process development.^13,14 Despite its importance a systematic understanding of confinement effects is still lacking.

Following up a recent work, in which we have investigated the thermodynamics and mechanism of phase transition between CO₂ forms I and III in bulk, we set out to study the effect of confinement on CO₂ condensed phases. We tackle this problem by carrying out a systematic analysis of CO₂ phase behaviour confined in weakly interacting cylindrical nano pores.

Cylindrical confinement is known to markedly affect phase behaviour. Even in systems of spherical particles, both rigid and soft, cylindrical nano-confinement has been found to induce characteristic packing structures, as well as abrupt transitions from order to disorder as a function of the void fraction.^15–20

Our work has a two-fold aim: on the one hand understanding phase behaviour of confined CO₂ is relevant due to its prominent role within the carbon cycle, and the surging needs for mitigating its emissions in atmosphere by implementing capture and storage technologies based on adsorption in porous solids.^21–24 On the other hand due to its modest structural complexity accompanied with a rich phase diagram, CO₂ represents a convenient model system to perform extensive sampling of polymorphic transitions at finite temperature²⁵ and gain insight on general aspects of molecular phase transitions under confinement.

Both experiments and theory highlight remarkable effects of confinements on CO₂ phase behaviour. For instance, the density of confined CO₂ can significantly exceed the fluid phase bulk density^26–31 reaching values comparable to that of solid CO₂ and displaying signs of orientational correlations and structural rearrangements.^26,27

Several examples of molecular modelling studies can be found investigating the transport of CO₂ in a confined fluid phase,^30,32,33 and the modelling literature characterising structural features of CO₂ under confinement is limited. An insightful work from Elola and Rodriguez has recently proposed a detailed analysis of the CO₂ fluid structure within silica pores, highlighting how confinement induces a layered arrangement in the fluid phase and determines a significant slowdown of both translational and rotational diffusion of CO₂ molecules.³¹

In our work we aim complementing the state of the art by understanding whether confinement promotes or inhibits the formation of ordered phases. To this aim we systematically analyse the intermolecular structural organisation of CO₂ confined in a simplified pore geometry consisting of a cylinder of Lennard-Jones (LJ) particles. We carry out a systematic study as a function of the cylinder radius (spanning from 1 to 5 nm), the nominal density of CO₂ molecules (5–15 molecules per nm³), and the interaction potential with the pore walls. We investigate this parameter space by systematically carrying out unbiased molecular dynamics simulations starting from liquid-like initial conditions, and by enhancing the exploration of ordered molecular packings with metadynamics.

Our analysis shows that, while the formation of known solid CO₂ arrangements is inhibited, confinement unveils a rather complex behaviour. Depending on the location in parameter space, confined CO₂ can either approach a completely disordered fluid state, exhibit short-range orientation correlations or arrange into long-range ordered packings. By systematically analysing relative orientations distributions we have systematically detected ordering events from molecular trajectories and classified ordered configurations.

This paper is organized as follows: firstly, in section 2, we introduce the model system, methods, parameters space, and the analysis tools developed to characterise ordering transitions. In section 3 we present and discuss our results, and finally in section 4 we summarise our findings by proposing a qualitative phase diagram in the space defined the nominal density of the pore and its normalised radius.

2 Methods

To investigate the effect of confinement on liquid carbon dioxide we have employed both standard molecular dynamics (MD)^34–36 and well-tempered metadynamics (WTMetaD).³⁷ For a detailed description of metadynamics we refer to Barducci et al.,^37,38 and Valsson et al.,³⁹ and for a brief overview of its applications in crystallisation studies to Giberti et al.⁴⁰ In the following we report the setup of the molecular model investigated, the details of MD simulations and WTmetaD simulations setup, and the analysis approach implemented to detect and characterise the emergence of an ordered phase from an initial liquid state.

CO₂ potential

In order to model CO₂ we employ the transferable potentials for phase equilibria (TraPPE) force field^41,42 with introduction of dummy atoms. As we aim at keeping the molecules rigid, this strategy, previously employed by Sanghi et al.³⁰ to represent 3-site CO₂, allows to avoid instabilities in MD simulations caused by constraining planar angles. The TraPPE force field is chosen among the models available for CO₂ for two reasons: on the one hand, it works reasonably well in the representation of the fluid at high pressure as well as of the liquid–solid border;^43,44 on the other hand, this choice is consistent with our previous work on carbon dioxide polymorphism at high pressure,²⁵ and allows for a direct comparison with results obtained in bulk.

Pore model

To develop a simplified model to study confinement in a weakly interacting porous medium, we build a single-wall cylindrical nanopore of LJ particles. We consider different pore diameters (1, 1.3, 2, 5 nm) with a constant wall density of 33.104 atoms per nm² and a height of 10 nm for all diameters except the case with diameter 5 nm, where the height has been limited to 5 nm for the sake of computational efficiency. The selection of pore diameters is based on literature works on confined CO₂;^27,31,45,46 we did not include larger sizes, as for values bigger than 2 nm we observe a bulk-like behaviour of the fluid. The cylinder axis aligns with z. We also investigate the effect induced by variations in the wall potential induced by keeping constant ε_wall (0.9977 kJ mol⁻¹) while varying the σ_wall, in the interval between 0.253 to 0.405 nm.

Simulation set up

For our simulations we choose temperature conditions of 323 K and nominal carbon dioxide densities compatible with supercritical CO₂.^24,28,30 The position of the wall atoms is frozen, but their intermolecular interaction with carbon dioxide molecules is taken into account, and, as presented in the results, plays a major role in the organization of the fluid inside the pore; CO₂ particles are instead allowed to move.

Firstly, we run NVT molecular dynamics simulations of the described system, with Bussi–Donadio–Parrinello⁴⁷ thermostat. We apply periodic boundary conditions (pbc) to ensure continuity along z, i.e. the cylinder axis. To avoid overlapping effects in the radial direction from periodic images, the simulation box is much larger (around 10 times) than the pore diameter. Lennard-Jones potential with Lorentz-Berthelot combination rules and long-range corrections is employed. Typical MD simulations are 20 ns-long, starting from disordered configurations equilibrated for 100 ps.

We use the same setup to carry out WTMetaD simulations aimed at exploring the space of accessible packings beyond the timescale limit of unbiased MD. As collective variable we employ an order parameter, hereafter indicated with λ, that expresses the degree of crystallinity of a system based on the local environment around each molecule and is able to distinguish not only between liquid and solid phases but also among polymorphs. A complete description of the mathematical formulation of λ has been reported by Giberti et al.⁴⁸ Explorative WTMetaD simulations are carried out biasing two distinct formulations of the λ order parameter: λ_I, tuned to capture the molecular arrangement of CO₂ bulk form I, which we have recently developed to investigate CO₂ polymorphism at high pressure;²⁵ and λ_B, based instead on the characteristic angles and local density of the most abundant ordered structure observed under confinement (configuration B in Table 2). Further explanation of the formulation and parameters for λ_I and λ_B, and details about the WTMetaD set up are available in the ESI.†

Parameter space representation

In order to represent our results throughout parameter space in a rational and efficient form we define a single adimensional radius, r′, which accounts for variations in both LJ potential (σ_wall) and pore size and quantifies the void space available within the pore.⁵ The adimensional radius r′ is defined by the following expression:where σ_m,C–wall represents the distance from the wall to the minimum of the C-wall potential (as in Fig. 1(a), for σ_wall = 0.34 nm) and σ_m,C–C is the position of the minimum of the C–C interaction. Values of r′ for all simulations performed in this work are reported in Table 1. The second parameter is the nominal density of CO₂ molecules within the pore, ρ_CO2, computed accounting for the volume of the cylinder of diameter d_pore.

Fig. 1 Simulation setup. (a) Representation of the LJ potential for the interaction between the pore wall and the carbon atom of CO₂; we report three pore sizes (d_pore = 1, 2, 5 nm), for which the zero on the x axis represents the cylinder axis and the walls are located at ±0.5, ±1 and ±2.5 nm, respectively; σ_wall is set to 0.34 nm. Arrows highlight σ_m,C–wall. (b) Intermolecular LJ potential characterising CC, CO and OO non bonded interactions. The x axis represents the interatomic distance; the black arrow indicates σ_m,C–C. (c) Representation of the cross-section of the initial configuration for simulations with d_pore = 1, 2, and 5 nm.

Table 1 Value of r′ for the combinations of d_pore and σ_wall investigated

r′		d pore [nm]
		1	1.3	2	5
σ wall [nm]	0.253	0.639	1.116	2.230	7.003
	0.315	0.528	1.006	2.119	6.892
	0.340	0.484	0.961	2.075	6.847
	0.372	0.472	0.904	2.017	6.790
	0.405	0.368	0.845	1.959	6.731

---

Systematic detection of ordered arrangements

In the definition of ordered parameters such as λ the identification of characteristic relative orientations is key.^25,40,49 This is typically based on the analysis of the distribution of relative orientations that characterises a relevant structure, as extensively discussed for bulk CO₂ in ref. 25. Indeed the relative orientation between molecular axes of neighbouring molecules represents a fingerprint of each arrangement and in principle allows to distinguish not only between liquid and solid, but also among different solid structures. An example of such distributions for bulk dry ice and liquid CO₂ is reported in Fig. 2. In the following we build on this observation to systematically detect ordering phenomena in within molecular trajectories.

Fig. 2 (a) The distribution of relative orientations is built considering all the angles between an assigned molecule and its nearest neighbours, falling within a sphere of radius r_cut represented in transparent green. The relative orientation is computed as the elevation angle between the axis of two CO₂ molecules. (b) Characteristic angle distribution for bulk phases, in particular melt (green) and phase I (blue); the analytical expression for the random distribution in melt is also reported in black.

As reported in Fig. 2, crystal and liquid carbon dioxide display significantly different characteristic orientation distributions. In order to quantitatively capture such difference we apply the Bhattacharyya distance,⁵⁰ a metric that quantifies the dissimilarity of two probability densities. The Bhattacharyya distance D_B between two angle distributions p(θ) and q(θ) is defined as:

D_B(p, q) = −ln(C_B(p, q)) (2)

where C_B is the Bhattacharyya coefficient, which measures the overlap between p(θ) and q(θ) defined as:

In our analysis we take advantage of the fact that a completely random arrangement of linear molecules displays a relative orientation probability density characterised by the functional form p_ref = 1/2 sin θ. For a derivation of the reference distribution functional form see the ESI.† As shown in Fig. 2 the distribution of angles in the bulk of liquid CO₂ closely resembles the theoretical random distribution p_ref.

Given a molecular trajectory, we can therefore define a time dependent measure of the deviation from an ideal random packing of CO₂ molecules as:

where p(θ,t) is the instantaneous probability density of relative orientations.

The construction of p(θ,t) is straightforward for carbon dioxide, due to its linear, symmetrical geometry. For each molecule we evaluate the polar angle between its molecular axis and the analogous vector in each one of its nearest neighbours, i.e. molecules within a cut-off radius of 4 Å for ordered phases and 5 Å for the liquid as shown in Fig. 2(a).

Tools

Pores are built using an in-house FORTRAN code, while we use Packmol⁵¹ to insert carbon dioxide molecules at random positions and orientations in order to avoid any creation of patterns. For MD and WTMetaD simulations we employ Gromacs 5.2.1⁵² and Plumed 2.2.⁵³ Post-processing of the results is carried on with MATLAB (R2015a), Visual Molecular Dynamics (VMD)⁵⁴ and Plumed 2.2.⁵³

3 Results and discussion

In this section we present and discuss the results of our molecular dynamics investigation of the effect of confinement on CO₂ ordering. Preliminary simulations were carried out to investigate the equilibrium distribution of a single CO₂ molecule in the presence of a cylindrical LJ wall and are discussed in ESI.† In the following we begin by discussing the arrangement of CO₂ molecules within the pore as a function of pore size and density, and we continue by reporting an analysis of the ordered arrangements identified within our parameter space.

Molecular arrangement within the pore

We begin by analysing the arrangement of CO₂ molecules in the pores at increasing ρ_CO2. In the case of low densities we observe the formation of a single layer of adsorbed CO₂ molecules in contact with the pore surface associated with a hollow structure. With increasing ρ_CO2 larger pores (d_pore = 2 and 5 nm), can accommodate a second layer of CO₂, and then display bulk-like filling as shown in Fig. 3(a and c). This finding is in agreement with the results obtained by Elola and Rodriguez³¹ for CO₂ confined in silica pores. Furthermore, the identified arrangement is reminiscent of the layered structure, identified by Prestipino et al.,²⁰ that emerges in a system of nanoconfined soft spheres due to the geometrical frustration of the freezing transition.

Fig. 3 Probability density profiles associated with the position of CO₂ molecules in the radial direction of the pore, for the range of densities investigated at different values of r′. In particular, (a) refers to r′ = 2.075, and (b) to r′ = 0.484. Along the abscissa, i.e. the pore radius, zero corresponds to the cylinder axis, while the maximum value reported to the radius of the pore. For the case in (b), (e) further analyses the height of the adsorbed layer peak and the number of peaks as a function of ρ_CO2. For the case in (a), instead, (c) reports snapshots of hollow to filled structures for increasing ρ_CO2, from 5.41 to 8.15 to 12.73 molecules per nm³. More in detail, (d) presents the probability density of the angle a as a function of the radial distance, for r′ = 2.075 and ρ_CO2 = 12.73 molecules per nm³. The angle α, reported as (cos α),² is evaluated between CO₂ molecular axis and the vector normal to the pore walls; to enhance the differences on the plot the probability density is represented by its logarithm. The x-axis convention is the same as before, with cylinder axis in zero.

Furthermore, we analyse CO₂ molecular orientation within the pore by evaluating the angle α between the molecular axis and the normal to the pore surface. As shown in Fig. 3(d), we can clearly note the high probability regions in the r, (cos α)² plane, corresponding to individual layers. We note that in both layers α is narrowly distributed around 90°. This indicates that CO₂ molecules arrange parallel to the pore surface, in agreement with observations reported in the literature.^27,30,31,55

In smaller pores (d_pore = 1 and 1.3 nm) a layered structure is not observed and the decrease in height of the dominant peak (peak 1 in Fig. 3(b)), accompanied by the appearance of a secondary peak (peak 2 in Fig. 3(b)) is instead associated with a reorganisation of the CO₂ packing within the pore (see Fig. 3(b and e)). For instance, we note that at d_pore = 1.3 nm, CO₂ arranges as a chain along the z axis for densities equal or above 9.87 molecules per nm³.

Detecting emerging order

After discussing the position and alignment of CO₂ molecules with respect to the pore wall, we delve into the analysis of intermolecular arrangement. In particular, we analyse the arrangement of the local environment around each CO₂ molecule by monitoring the distribution of relative orientations. As mentioned in section 2, this information constitutes a structural fingerprint and can discriminates not only between ordered and disordered arrangements but also among different ordered packings.

As discussed in section 2, to detect the unfolding of ordering events we compute the Bhattacharyya distance, D_B(t) with respect to an ideal random packing for each MD trajectory. When ordered arrangements nucleate within a pore, the D_B metric displays a sharp increase. An example of this behaviour in a typical trajectory is reported in Fig. 4, where in (b) we compare the trend of D_B(t) for a trajectory undergoing an ordering transition (violet) and a trajectory that remains stable in the fluid state (green). In Fig. 4(a) the evolution of the instantaneous intermolecular angle distribution p(θ,t) is reported together with the p_ref(θ) corresponding to an ideal random packing.

Fig. 4 Bhattacharyya analysis of an ordering phenomenon. (a) Temporal evolution of the angle distribution, p(θ,t), for r′ = 0.484 − ρ_CO2 = 11.78 molecules per nm³, where colour dark blue to light blue represents time increment, while the black distribution is the probability density for melt, employed as reference p_ref for evaluating D_B. (b) Bhattacharyya distance, D_B, as a function of time for an example of ordering (violet, r′-ρ_CO2 as (a)) and fluid (green, r′ = 0.961 − ρ_CO2 = 5.41 molecules per nm³ conditions) the probability density of D_B is also reported in the form of a histogram.

Mapping phase behaviour in the (ρ_CO2, r′) parameter space

The systematic analysis of the deviation from an ideal random packing and the contextual detection of the appearance of ordered phases from a liquid configuration allows to clearly map the qualitative behaviour of confined CO₂ within the (ρ_CO2, r′) parameter space. In particular we can identify areas where CO₂ under confinement spontaneously evolves to ordered packings within the characteristic timescale of MD simulations. To further analyse our finding, we complement the analysis of the angle distribution by reporting the C–C radial pair distribution function g(r). Inspecting both g(r) and relative orientation distribution we can identify three regions in parameter space: a disordered region (Fig. 5 green area, examples (a and b)), an ordered region (Fig. 5 violet area, examples (e–f)), and a transition region between them (Fig. 5 green-to-violet area, examples (c and d)). We note that in the disordered region, at large r′ and low density, the distribution of relative orientations approaches the ideal random distribution and the g(r) displays the short-range order typical of a liquid. On the contrary, in the ordered region both the angle distribution and the g(r), deviate from the typical liquid behaviour denoting long-range order through the presence of characteristic peaks in both orientation distribution and g(r). Finally, the transition region denotes a peculiar behaviour. In this region of parameter space, while g(r) maintains a close resemblance to that of a liquid phase, the distribution of relative orientations displays significant deviations from an ideal random packing suggesting the development of correlation in the relative orientation of CO₂ molecules.

Fig. 5 Radial pair distribution function, g(r) (black) and angle distribution (blue) for illustrative cases of liquid (a and b), transition (c and d), and ordered structures (e and f). In particular, (a) refers to r′ = 2.075 − ρ_CO2 = 5.41 molecules per nm³, (b) to r′ = 2.075 − ρ_CO2 = 12.73 molecules per nm³, (c) to r′ = 0.484 − ρ_CO2 = 5.41 molecules per nm³, (d) to r′ = 0.484 − ρ_CO2 = 8.15 molecules per nm³, (e) to r′ = 0.484 − ρ_CO2 = 11.78 molecules per nm³, and (f) to r′ = 0.484 − ρ_CO2 = 12.73 molecules per nm³. Green areas on the phase diagram refer to fluid state, violet to ordered, while violet-to-green colour to transition structures.

From the analysis of the transitions obtained in the structured region of parameter space we notice a finer grained level of structural complexity. As displayed by the qualitative comparison of both the distribution of relative orientations and the g(r) reported in examples Fig. 5(e) and (f) multiple ordered packings can be obtained.

Stable ordered packings

We now move on to present a detailed analysis of the ordered structures encountered within the ordered region of the (ρ_CO2, r′) phase diagram. We shall note that, in order to confirm the finding emerging from unbiased simulations, we have applied WTMetaD simulations to further sample the configuration space of confined CO₂. This allows to verify that in the disordered or transition regions ordered arrangements cannot be found or, when found do not states with a finite lifetime and spontaneously revert to liquid in unbiased MD.

Nevertheless, in the region of the phase diagram where confinement-induced order emerges, we identify four different stable arrangements. We assign to such arrangements a label progressing from A to D, based on their structural features. As we can see in Table 2, progressing from A to D these structures present an increasing number of molecules in their repeating unit. Moreover, the same classification applies when analysing the geometry of the apparent ring that characterises the horizontal cross-section of the pore (see Table 2). The number of ring members is 3 in A, 4 in B, and 6 in both C and D, with D displaying one molecule in the middle of the pore oriented parallel to the pore longitudinal axis.

Table 2 Molecular representation of the arrangement and characteristic angles distribution for the ordered structures identified. Such structures are named A to D according to the apparent number of ring members when observed from the top. The angles distributions are obtained at Γ = 0.484 − ρ_CO2 = 9.87, at r′ = 0.484 − ρ_CO2 = 11.78, at r′ = 0.484 − ρ_CO2 = 12.73, at r′ = 0.961 − ρ_CO2 = 12.73, for A, B, C, and D, respectively

Label	Snapshot	Relative orientation distribution
A
B
C
D

---

As far as the packing of these four ordered structures is concerned, it can be seen that the characteristic angle distribution represents a clear fingerprint of their molecular arrangement (Table 2), while the coordination number appears to be a function of the total density inside the pore and cannot clearly discriminate between different arrangements. In particular, the orientation distribution of B presents well-defined peaks in θ₁ = 46.2° and symmetrically in θ₂ = 133.8°, while C a very sharp single one in θ = 90°. Structure A, instead, appears to be an intermediate between B and C both from the inspection of the coordinates and from its characteristic angle distribution. Indeed, while the most densely populated orientations closely resemble those of B, the peak at 90° is more pronounced than in B. Form D present instead wider peaks in θ₁ = 57.9° and the additional peak θ₂ = 122.1°.

Interestingly, three out of four ordered structures, more precisely A, B, and D, form spontaneously in unbiased MD simulations; phase C, instead, is generated from metadynamics simulations performed in conditions that spontaneously yielded other ordered arrangements (violet area in the r′/ρ_CO2 phase diagram in Fig. 5 and 6); structure C is confirmed to be a metastable state for the system possessing a finite lifetime through 20 ns long unbiased MD simulations, initialised in configuration C, which do not display any sign of instability transition.

It should be noted that the exploration carried out with WTMetaD yields additional ordered arrangements, which, however, do not survive unbiased MD and are therefored considered to be unstable (see ESI†).

Fig. 6 Qualitative phase diagram of CO₂ confined in weakly attractive cylindrical nanopores. In green is represented the region of phase space characterised by disordered configurations, while we report in violet the region of phase space where ordered structures are found. Configurations yielding one or more amongst packing A, B, C, and D are appropriately marked. Qualitative phase boundaries have been reported in order to highlight how simulations yielding similar structures tend to cluster in parameter space, suggesting that the spontaneous ordering process observed from MD reflects the underlying thermodynamic stability of the packings.

In addition metadynamics further confirms that for ρ_CO2-r′ conditions that induce spontaneous ordering, the disordered liquid state is strongly disfavoured with respect to the ordered phases. Indeed even in extended WTMetaD simulations, characterised by the deposition of significant bias, no melting event can be observed, while transitions between ordered arrangements are routinely detected. This finding is consistent irrespective of the set of CVs used to enhance the exploration and suggests that the energetic barrier associated with melting is higher than the barriers associated with solid–solid interconversion.

Furthermore, the conversion from hollow (A, B, C) to filled (D) ordered structures and vice versa is not sampled during WTMetaD simulations; this transition appears unlikely in the conditions investigated since, as discussed before, different values of r′ can accommodate different numbers of CO₂ layers.

On the other hand, transition between phase A, B and C is possible and consistently sampled. More precisely, it takes place though a gradual reorganisation of the layers that starts in a localised region and then spreads through the entire pore length. Interestingly, WTMetaD allows to sample mixed arrangements of different packings, (e.g. half ordered as B and half as C) which reveal to be stable in subsequent MD simulation, confirming that solid–solid interconversion is an activated event. Moreover, in these structures it is common to find boundary regions in which the order is not well-defined, effectively acting as defects within the pore.

4 Conclusions

In this work we present a systematic investigation of the effects of confinement in cylindrical nanopores on CO₂, as a function of pore size, wall potential and CO₂ density. To systematically detect the formation of ordered packings we have devised a general analysis approach based on the calculation of the time-dependent Bhattacharyya distance between a reference probability density and the instantaneous probability density of the relative orientation of nearest neighbours. Here we have adopted as a reference an ideal random packing, which closely approximates the distribution of relative orientations in the bulk of the liquid phase. This analysis approach is general and we anticipate its applicability to a much wider range of phase transitions in molecular solids.

We find that cylindrical confinement induces polymorph selection by inhibiting the nucleation of known solid forms of CO₂ while at the same time inducing the organisation of CO₂ into a series of distinct ordered arrangements, which do not resemble any of the known CO₂ polymorphs. Using a combination of unbiased MD and WTMetaD we identify four ordered packing that are stable in the region of parameter space corresponding to small pore radii and large densities. These configurations have been labelled A, B, C and D following the progression in the number of CO₂ molecules visible in a cross-section of a pore. Interestingly, none of these configurations corresponds to known bulk structures. It should be noted that, while A, B, and D emerge spontaneously from unbiased MD, C is found from the enhanced exploration of the configuration space achieved with WTmetaD.

In Fig. 6 we report a qualitative phase diagram that summarizes the identified structures and the regions of parameter space where such structures have obtained. Markers point the location in parameter space where specific packings have been observed. Cases in which multiple configurations are indicated on the phase diagram represent points in which multiple arrangements have been observed either with unbiased MD or through WTMetaD. We observe that, in the ordered region of the phase diagram WTMetaD was quite efficient in promoting transitions between ordered packings, however it was unable to induce melting. Despite its qualitative character, this observation suggests that within the ordered region of the phase diagram solid–solid interconversion is energetically favoured with respect to liquid mediated transitions.

Finally we observe that, even in the disordered region confinement induces a short-range organisation of the liquid phase, in which the distribution of relative orientations departs from an ideal random packing, as shown in Fig. 5. At large volumes and low densities the bulk behaviour is recovered as expected.

To conclude we note how such a simple model system such as CO₂ confined in a simple nanometric pore model unveils remarkable structural complexity, showcasing the profound effect of confinement on the phase behaviour of molecular solids. This work paves the way for a systematic assessment of CO₂ packing polymorphism in systems characterised by realistic pore geometries and surface chemistry.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge EPSRC (Engineering and Physical Sciences Research Council) for PhD scholarship (EPSRC Award ref 1760926), and UCL Legion High Performance Computing Facility for access to Legion@UCL and associated support services, in the completion of this work.

References

1. J.-M. Ha, J. H. Wolf, M. A. Hillmyer and M. D. Ward, J. Am. Chem. Soc., 2004, 126, 3382–3383.
2. A. Llinàs and J. M. Goodman, Drug Discovery Today, 2008, 13, 198–210.
3. Y. Maniwa, H. Kataura, M. Abe, S. Suzuki, Y. Achiba, H. Kira and K. Matsuda, J. Phys. Soc. Jpn., 2002, 71, 2863–2866.
4. A. Striolo, A. A. Chialvo, K. E. Gubbins and P. T. Cummings, J. Chem. Phys., 2005, 122, 234712.
5. H. Kyakuno, K. Matsuda, H. Yahiro, Y. Inami, T. Fukuoka, Y. Miyata, K. Yanagi, Y. Maniwa, H. Kataura, T. Saito, M. Yumura and S. Iijima, J. Chem. Phys., 2011, 134, 244501.
6. G. Algara-Siller, O. Lehtinen, F. C. Wang, R. R. Nair, U. Kaiser, H. A. Wu, A. K. Geim and I. V. Grigorieva, Nature, 2015, 519, 443–445.
7. A. K. Soper, Nature, 2015, 519, 417–418.
8. K. V. Agrawal, S. Shimizu, L. W. Drahushuk, D. Kilcoyne and M. S. Strano, Nat. Nanotechnol., 2016, 12, 267–273.
9. J. Chen, G. Schusteritsch, C. J. Pickard, C. G. Salzmann and A. Michaelides, Phys. Rev. Lett., 2016, 116, 025501.
10. J. Chen, A. Zen, J. G. Brandenburg, D. Alfè and A. Michaelides, Phys. Rev. B, 2016, 94, 220102.
11. C. J. Stephens, S. F. Ladden, F. C. Meldrum and H. K. Christenson, Adv. Funct. Mater., 2010, 20, 2108–2115.
12. Y.-W. Wang, H. K. Christenson and F. C. Meldrum, Adv. Funct. Mater., 2013, 23, 5615–5623.
13. A. Y. Lee, I. S. Lee, S. S. Dette, J. Boerner and A. S. Myerson, J. Am. Chem. Soc., 2005, 127, 14982–14983.
14. M. Beiner, G. T. Rengarajan, S. Pankaj, D. Enke and M. Steinhart, Nano Lett., 2007, 7, 1381–1385.
15. R. O. Erickson, Science, 1973, 181, 705–716.
16. F. Durán-Olivencia and M. Gordillo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 061111.
17. M. Maddox and K. Gubbins, J. Chem. Phys., 1997, 107, 9659–9667.
18. M. A. Lohr, A. M. Alsayed, B. G. Chen, Z. Zhang, R. D. Kamien and A. G. Yodh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 040401.
19. A. Mughal, H. Chan, D. Weaire and S. Hutzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051305.
20. S. Prestipino, F. Saija, A. Sergi and P. V. Giaquinta, Soft Matter, 2013, 9, 9876–9886.
21. S. M. Benson, in Carbon Dioxide Capture for Storage in Deep Geologic Formations, Elsevier, 2005, vol. 2, pp. 665–672.
22. F. M. Orr, Science, 2009, 325, 1656–1658.
23. N. MacDowell, N. Florin, A. Buchard, J. Hallett, A. Galindo, G. Jackson, C. S. Adjiman, C. K. Williams, N. Shah and P. Fennell, Energy Environ. Sci., 2010, 3, 1645–1669.
24. S. Krevor, M. J. Blunt, S. M. Benson, C. H. Pentland, C. Reynolds, A. Al-Menhali and B. Niu, Int. J. Greenhouse Gas Control, 2015, 40, 221–237.
25. I. Gimondi and M. Salvalaglio, J. Chem. Phys., 2017, 147(11), 114502.
26. G. K. Papadopoulos, J. Chem. Phys., 2001, 114, 8139–8144.
27. T. Steriotis, K. Stefanopoulos, N. Kanellopoulos, A. Mitropoulos and A. Hoser, Colloids Surf., A, 2004, 241, 239–244.
28. Y. B. Melnichenko, H. Mayama, G. Cheng and T. Blach, Langmuir, 2010, 26, 6374–6379.
29. G. Rother, E. G. Krukowski, D. Wallacher, N. Grimm, R. J. Bodnar and D. R. Cole, J. Phys. Chem. C, 2012, 116, 917–922.
30. T. Sanghi and N. R. Aluru, J. Chem. Phys., 2012, 136, 024102.
31. M. D. Elola and J. Rodriguez, J. Phys. Chem. C, 2016, 120, 1262–1269.
32. X. Yang and C. Zhang, Chem. Phys. Lett., 2005, 407, 427–432.
33. A. Botan, B. Rotenberg, V. Marry, P. Turq and B. Noetinger, J. Phys. Chem. C, 2010, 114, 14962–14969.
34. M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, OXFORD UNIVERSITY PRESS, 2010.
35. D. Frenkel and B. Smit, Understanding Molecular Simulation-From Algorithms to Applications, ACADEMIC PRESS, 2nd edn., 2002.
36. J. C. Palmer and P. G. Debenedetti, AIChE J., 2015, 61, 370–383.
37. A. Barducci, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 1–4.
38. A. Barducci, M. Bonomi and M. Parrinello, WIREs Comput. Mol. Sci., 2011, 1, 826–843.
39. O. Valsson, P. Tiwary and M. Parrinello, Annu. Rev. Phys. Chem., 2016, 67, 159–184.
40. F. Giberti, M. Salvalaglio and M. Parrinello, IUCrJ, 2015, 2, 256–266.
41. J. J. Potoff, J. R. Errington and A. Z. Panagiotopoulos, Mol. Phys., 1999, 97, 1073–1083.
42. J. J. Potoff and J. I. Siepmann, AIChE J., 2001, 47, 1676–1682.
43. C. G. Aimoli, E. J. Maginn and C. R. Abreu, Fluid Phase Equilib., 2014, 368, 80–90.
44. G. Pérez-Sanchéz, D. Gonzaléz-Salgado, M. M. Piñeiro and C. Vega, J. Chem. Phys., 2013, 138, 084506.
45. E. Pantatosaki, D. Psomadopoulos, T. Steriotis, A. Stubos, A. Papaioannou and G. Papadopoulos, Colloids Surf., A, 2004, 241, 127–135.
46. S. K. Bhatia, K. Tran, T. X. Nguyen and D. Nicholson, Langmuir, 2004, 20, 9612–9620.
47. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
48. F. Giberti, M. Salvalaglio, M. Mazzotti and M. Parrinello, Chem. Eng. Sci., 2015, 121, 51–59.
49. M. Salvalaglio, T. Vetter, F. Giberti, M. Mazzotti and M. Parrinello, J. Am. Chem. Soc., 2012, 134, 17221–17233.
50. A. Bhattacharyya, Bull. Calcutta Math. Soc., 1943, 35, 99–109.
51. D. G. Truhlar, J. Comput. Chem., 2009, 28, 73–86.
52. M. Abraham, D. van der Spoel, E. Lindahl, B. Hess and the GROMACS development Team, GROMACS User Manual version 5.1.1, www.gromacs.org, 2015.
53. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
54. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
55. A. Vishnyakov, P. I. Ravikovitch and A. V. Neimark, Langmuir, 1999, 15, 8736–8742.

---

Footnote

† Electronic supplementary information (ESI) available: Single molecule probability distribution in the model pore, derivation of the ideal random distribution p_ref(θ), details of the setup and results of explorative metadynamics simulations, additional analysis of the radial density profiles. See DOI: 10.1039/c7me00103g

---