Atomistic Monte Carlo Simulations on the Formation of Carbonaceous Mesophase in Large Ensembles of Polyaromatic Hydrocarbons

Nano-carbons comprise a class of advanced materials generally associated with high strength and light weight with applications as composites, heat-sinks and battery electrodes, and in the aerospace, automotive and sports-equipment industries [1]. On the molecular scale, carbon materials are composed largely of graphene layers – planar layers of carbon atoms which are covalently bonded into hexagonal aromatic rings [2]. These graphene layers pack together and fold over one another to form a variety of different microscopic textures. The amount of order present in this texture greatly influences the macroscopic mechanical properties of the material, and so is an important consideration for the materials engineer. A well-ordered texture, for example, can be readily converted into an anisotropic, graphitic structure (in a process known as graphitisation), whereas a highly disordered texture tends to simply fuse into an isotropic, glassy structure [3]. The clear distinction between these two types of carbons (graphitising and non-graphitising respectively) was originally made by Frankin in 1951 [4] based on the X-ray diffraction of several carbon materials.


Introduction
Nano-carbons comprise a class of advanced materials generally associated with high strength and light weight with applications as composites, heat-sinks and battery electrodes, and in the aerospace, automotive and sports-equipment industries [1]. On the molecular scale, carbon materials are composed largely of graphene layers -planar layers of carbon atoms which are covalently bonded into hexagonal aromatic rings [2]. These graphene layers pack together and fold over one another to form a variety of different microscopic textures. The amount of order present in this texture greatly influences the macroscopic mechanical properties of the material, and so is an important consideration for the materials engineer. A well-ordered texture, for example, can be readily converted into an anisotropic, graphitic structure (in a process known as graphitisation), whereas a highly disordered texture tends to simply fuse into an isotropic, glassy structure [3]. The clear distinction between these two types of carbons (graphitising and non-graphitising respectively) was originally made by Frankin in 1951 [4] based on the X-ray diffraction of several carbon materials.
Carbon materials are often derived from the poly aromatic hydrocarbons, a major component of pitch materials. 80% of industrial carbon is produced from pitch precursors, particularly from petroleum-derived pitches [5]. Pitch is a cheap material, and its low cost gives it a wide range of potential applications. When a pitch is converted into a carbon material, the degree of textural ordering of the graphene layers depends strongly on the formation of a mesophase at high temperatures. A range of carbonaceous materials are formed from pitch precursors, the nanostructures of which are dependent upon the formation of this mesophase. Examples of these include carbon fibres, nanofibres, nanotubes, and foams.
One of the key issues during the processing of carbon materials is the formation carbonaceous mesophase, a liquid crystalline phase that can form during the heat treatment of a variety of pitch materials [6,7]. As the anisotropy of the final product is determined to some extent by this phase, the formation or otherwise of this phase is a key step in the production of longrange ordered carbon structures. It therefore becomes important to understand the role played by factors such as molecular size, shape, the presence of aliphatic rings etc. driving this phase transformation [8].
Despite the importance of this phase, little is known about its fundamental atomistic behaviour. It is currently very difficult to obtain polycyclic aromatic hydrocarbons in sufficient purity for experimental studies of the phase behaviour [9]. Pitches generally contain a large distribution of molecule sizes/shapes and the isolation of a specific molecule is very difficult unless these molecules are artificially grown under controlled conditions [10]. Another key difficulty is that the polymerisation of the molecules is unavoidable at the temperatures required to experimentally observe the mesophase, thereby destroying the original molecule being studied [11]. To the best of our knowledge, a direct fundamental experimental investigation on carbonaceous meso-phase is currently not possible on a molecular level.
Polycyclic aromatic hydrocarbons (PAH) of high molecular weight (> 400 Daltons) are the primary constituents of pitches [12]. These are formed in situ by cyclization, dehydrogenation, and polymerization during carbonization of diverse organic materials [13]. Although carbonaceous meso-phase has been experimentally observed in a number of systems, these are generally complex mixtures of thousands of molecules. The fundamental information regarding specific molecular constituents is difficult to come by from experimental data. Molecular weights of about 600 or more are believed to be required for the formation of carbonaceous mesophase, the discotic nematic liquid crystalline intermediate that provides an effective route to long-range (> 1 μm) crystalline order in carbonaceous materials [14].
Molecular modelling and computer simulation can play an especially important role in such a study. There have been a number of molecular modelling studies focused on non-covalent PAH interactions [15][16][17], but the size and complexity of the molecules involved have restricted most studies to energy calculations on polyaromatic clusters (dimers, trimers, and tetramers) and simulations of large disk ensembles using molecular pseudo-potentials [18,19]. These large PAH molecules assemble by non-covalent interactions into supra-molecular structures. There have been a few interesting studies of Gay-Berne discotic phase behaviour [20], but there is no obvious way to relate detailed molecular structure (aliphatic groups, non-planar heterocylic rings, irregular polyaromatic cores) to the empirical potential parameters.
Our group has developed an isobaric-isothermal Monte Carlo technique for simulating large ensembles of discotic molecules based on massive summation of atomic pair potentials [21]. This technique explicitly uses description of atom locations with no adjustable potential parameters. A parametric study of molecular structural features such as non-planar aliphatic rings, irregularly-shaped polyaromatic cores was carried out. The formation of carbonaceous meso-phase was reported in an irregular pitch molecule containing 58C and 28H atoms. The phase behaviour was found to be sensitive to minor perturbations in the molecular structure.
A study was also carried out on a large regular molecule that also showed the formation of meso-phase [22].
In this paper, we report computer simulations results on the phase behaviour of five organic molecules in a range of sizes and shapes. The aim is to determine the effect of molecular size and shape on the formation of carbonaceous meso-phase towards developing an understanding of basic mechanisms and factors affecting this transition. The article is organised as follows.
In section 2, we present details of molecules under investigation and their specific features. Simulation algorithms, macromolecular interaction potentials and results on ordered ground state structures are presented in Section 3. Computer simulations at elevated temperatures for determining the phase behaviour will be presented in Section 4. A brief discussion on the mechanisms of phase transition and conclusions will be presented in Section 5.

Molecular structures
A number of criteria were used in the choice of model molecules.

Shape:
Each molecule was to have a shape which resembled a disc. Hexagonal rotational symmetry was enforced, keeping the molecules as close to a circular shape as possible. This shape contrasts with the irregular shapes of pitch mesogens studied previously with this technique allowing a quantitative comparison to be made based on molecular regularity. On the other hand, it closely resembles with parametric discotic mesogens studied in non-atomistic simulations.

Molecular weight:
The molecules were selected to represent a variety of molecular weights covering the range 400-1000 Da, the range typical for mesophase-forming pitches [23]. Several molecular weights lying close to and yet outside of this range were also selected to study critical limits of molecular weights.

Planarity:
It has been shown [21] that changes in molecular planarity can have a significant impact on mesophase formation. Aliphatic rings are not entirely planar, and necessarily introduce additional hydrogen atoms that are out of plane from the rest of the molecule. A fully aromatic structure, on the other hand, enforces planarity. One study was also carried out on an irregular shaped non-planar molecule containing an aliphatic ring.
Five molecules were investigated in this study. Their specific details are given below: a. Coronene: This regular shaped molecule contains seven aromatic rings containing 24C and 12H atoms on the edges of outer rings. Its molecular weight is 300 Daltons. It is one of the largest molecules that have been isolated and investigated experimentally.
b. Hexa-benzo coronene: This regular shaped molecule contains thirteen aromatic rings containing 42C and 18H atoms on the edges of outer rings. Its molecular weight is 523 Daltons. It can be constructed by adding six benzene rings on the outer edges of a coronene molecule.
c. Circum-coronene: This regular shaped molecule contains nineteen aromatic rings containing 54C and 18H atoms on the edges of outer rings. Its molecular weight is 667 Daltons. It can be constructed by adding 12 benzene rings on the outer circumference of a coronene molecule.
d. Dodeca-benzo-circumcoronene: This regular shaped molecule contains thirty one aromatic rings containing 84 C and 24H atoms on the edges of outer rings. Its molecular weight is 1033 Daltons. It can be constructed by adding 12 benzene rings on the outer circumference of a circum coronene molecule.
e. Representative pitch molecule: This irregular shaped molecule contains one aliphatic ring and seventeen aromatic rings. It contains 58C and 24H atoms. Its molecular weight is 720 Daltons. Some studies on this molecule were performed treating aliphatic ring as an aromatic ring to make the molecule planar in order to compare its behaviour with other four molecules under investigation.

Interaction potentials
A first step in investigating the high temperature phase behaviour of these molecules was to determine their ground state ordered structures as it is generally much easier to study orderdisorder transitions from the ordered state than creating order from disorder, especially in the absence of experimental data. These simulations were carried out in two steps. Assuming the molecule to be rigid in nature, the energetics of molecular packing into one and two-dimensional aggregates was described by non-bonded atom-atom potentials (E nb ) and an electrostatic contribution (E el ) from the partial charges on various atoms.
We have used the MM2 force field of Allinger [24] for these computations and atomistic modelling of macromolecules. A more recent version of these force fields, MM3, incorporating a slightly larger size of atoms and somewhat softer potentials as compared to MM2, was found to be too soft, with rather weak repulsion at small distances and not very suitable in these simulations. [25] The non-bonded potential for the interaction of two molecules labelled as 1 and 2 is given by 6 12 , 12.5 2.90 10 exp 2.25 where r ij is the distance between the atom i in the molecule 1 and the atom j in the molecule 2.
A ij and B ij respectively represent non-bonded van der Waals energy and distance parameters. An average C-C bond length of 1.41Å and an average C-H bond length of 1.09 Å was used in these simulations. For the C-H bond, the bond length r ij has to be shrunk by 8.5% before energy computations to account for the anisotropy of non-bonded interactions which depend on the relative orientations of interacting bond orbitals [26]. Because of π electron interactions, aromatic and saturated carbons need to be treated differently and have slightly different intermolecular potentials. The van der Waals interaction parameters for C and H atoms, based on MM2 semi-empirical potentials are listed in Table 1.  [20] The contributions from the electrostatic term, E el , to the total packing energy of molecules were neglected for these hydrocarbons that lack dipoles. The neglect of electrostatic terms here is not believed to be a major limitation in the technique, and indeed they could be added within this framework if necessary, such as for systems involving hetero-atoms.

2D Simulation of ordered structures
The ability of organic molecules to form an ordered crystalline array is one of the most remarkable occurrences in nature. While the crystal structure of a large number of organic solids has been measured experimentally, there is a great interest in predicting the structure of crystalline solids from first principles. This assumes a great significance for high-MW PAH, where a crystallographic investigation is not feasible for a number of reasons. Using the methodology developed by Perlstein [27], we have computed in two steps, a 2D structural arrangement for all molecules under investigation. In the first step, a lowest energy 1D aggregate structure was generated using Monte Carlo simulations. A 2D monolayer was then determined in the second step, consistent with lowest energy and closest packing. A Monte Carlo cooling technique was used for computing the global energy minimum and conformation of the molecule in different 1D aggregates [28].
One of the most remarkable features in organic crystallography is that only a very small number of 1D chains and 2D layers that actually occur in nature. Scaringe and Perez [29] have shown that only four types of 1D aggregates occured in 92% of all crystals: these include translation, screw (2 1 Screw), glide and inversion. Scaringe has shown that these can combine to form only seven types of layers. 1D aggregate containing five PAH molecules was constructed in four most probable symmetry arrangements. The binding energy of the aggregate corresponds to half the binding energy of the central molecule in the aggregate. In the initial configuration, the centre of mass of the molecule was placed at the origin of an orthogonal coordinate system and then all atoms in the molecule were put in place. This molecule, in an arbitrary orientation, was treated as the original. The z-axis was chosen as the axis of symmetry.
The energy of the aggregate was computed using equations (1) and (2). Using a random number generator, random changes to three orientation angles and translation/offset distances were generated within preset limits. The maximum range of these variables was determined by ± θ max , z min -z max and x min -x max . While θ max was chosen as 10°, the changes in the x and z distances were kept to a maximum of 6Å. A new aggregate was constructed using modified dimensions and the energy change was computed. The change in dimensions was accepted for ΔE ≤ 0. For ΔE > 0, the change could be accepted with a transition probability W, [30] defined as where k B and T are respectively the Boltzmann constant and temperature. A dimensionless, reduced temperature parameter, T* ( = k B T / A c-c ), was used to represent the simulation cell temperature. A c-c ( = 3.8 kcal/mol) represents the nonbonding C-C interaction parameter A from MM2 force fields (Table 1). W was compared to a random number, chosen uniformly between 0 and 1. The change was accepted for and the Monte Carlo step was repeated with the new configuration. The following cooling schedule was used in the search for the global minimum.
Using large values of angular and positional changes, Monte Carlo simulations were started at high temperature, e.g., T* = 4 and 50 cycles were carried out. The system was cooled to 0.9T* with the dimensional changes also reduced by a factor of 0.9 after which another 50 cycles were carried out. System cooling was continued in small steps until the temperature reached 0.1T. Treating the final state thus obtained as the new initial state, the simulations were restarted from step (1). The energy changes and downhill energy movements were monitored throughout. The entire process was repeated at least twenty times. Simulations were continued till no further decrease in energy could be observed. These simulations typically required approximately 1 hour of CPU time.
With continuous thermal cycling and downward energy movement, the system was expected to reach the lowest energy configuration at the end of simulations. While the final configuration was expected to be the global minimum, it is quite probable that the system may pass through other very low energy states during these simulation cycles. Whenever the energy of an aggregate went below a certain preset energy level, its molecular configuration and energy were recorded. It was then possible to establish the lowest energy configuration unambiguously in all cases. It can be seen that for coronene, the simulated Herring Bone structure with screw 1D symmetry compared very well with the experimentally determined structure for this molecule. This structure was also favoured by hexa-benzo coronene as well as by circum-coronene. Other two molecules favoured translational symmetry. The lowest energy configuration in 2D needs to satisfy the close packing criterion as well, these vertical columns were tilted with respect to the horizontal in order to conserve space while avoiding an overlap between the atoms of neighbouring molecules.

Simulation algorithm
An isothermal-isobaric Monte Carlo simulation approach was used to study the phases exhibited by the model PAH molecules as a function of temperature. The simulation lattice was a 16 x 16 monolayer, with the molecules occupying the ordered lattice in their ground state configuration. The outer ring of molecules in the simulation lattice was used as fixed walls, giving the cell a fixed volume. While the molecules in the outer ring were kept fixed and did not undergo any change in orientation or position during a simulation run, these participated in the energy computations of inner cell molecules. The application of pressure after a number of constant volume runs, allowed for the volume and shape changes to the simulation cell permitting the structure to attain its natural state. Pressure also contributed to system energy through a PΔV term, where P is the pressure parameter and ΔV is the change in the simulation cell volume. The details of the Monte Carlo algorithm are as follows.
A molecule was chosen at random from the inner 14 x 14 simulation lattice. Keeping the molecule rigid, it was allowed one of the five possible moves: three changes to the molecular orientation and two positional movements. Maximum ranges of angular rotation were: Δθ max : ± 3°, and displacement: horizontal motion ΔX max : ±1Å, and vertical motion ΔZ max : ±0.5Å. A number set between 0 and 1 was divided into five equal subsets, which were then allocated to specific molecular movements. A random number η was chosen uniformly between 0 and 1 and directed to the appropriate subset. The deviation of η from the midpoint of the subset, determined the sign and the magnitude of the specific movement. Simulations were also carried out by simultaneously applying multiple movements to a given molecule. These however showed a poor convergence with most of the moves becoming energetically unfavourable and were rejected.
The energy of the system was computed before and after the molecular movement. The computation of energies was the most time consuming part of the simulation algorithm. The cut-off distance for atomic interactions was kept at 15Å. Energy computations were carried out on the selected molecule within a 5x5 neighbour cell around it. A further increase in the size of the interaction cell was found quite unnecessary at all temperatures. The energy change due to the molecular move was calculated. The molecular movement was accepted for ΔE ≤ 0. For ΔE > 0, the move could be accepted with the Boltzmann transition probability.
The simulations were carried with the reduced temperature parameter T ranging from 0.1 to 1.5 and the pressure parameter P ranging from 1 to 10. Temperature and pressure were kept constant during a given run. Typically 100,000 to 500,000 Monte Carlo cycles were found sufficient ensuring the system had reached equilibrium. Each cycle consisted on ~2000 molecule trial moves at constant volume followed by a trial shape change, in which the dimensions of the simulation cell were randomly altered. Some of the simulations were also carried out where the final temperature was reached in two to three steps instead of a single jump from 0 to the final temperature. The final equilibrium configuration was analysed in terms of system energy, centre of mass/atomic positions and molecular orientations, and the orientational order parameter S (= 3 2 < cos 2 (θ) > − 1 2 ; where the brackets denote the average over all molecules). A number of pressure parameters were trialled, a parameter of 5.0 was found to be quite adequate for most simulations. These simulations can take up to 60 hours of CPU time. Most of the simulations used batch processing as these Monte Carlo algorithms cannot be efficiently transformed for parallel processing. Up to five alternatives were used near phase transformations and up to three alternatives were used for temperatures away from the transition. The phase transitions were approached from the high temperature and the low temperature sides and simulation time was much extended (up to 3 times) to account for the critical slowing down near a phase transition and an increased level of fluctuations.

Simulation results:
Simulation results from five model molecules are presented below.

Coronene
Carbonaceous mesophase was not observed in the case of coronene. This result was expected because of the small size of the molecule. It showed a complete transformation from the ordered phase to random orientational and positional disorder with increasing temperatures and simulation time. A typical energy plot of the system is shown in Figure 2.

Hexa-benzo Coronene
Hexa-benzo Coronene also showed a behaviour quite similar to that of coronene with the system going directly to complete disorder from the initially ordered state. The orientational order parameter also showed a continuous decrease from 1 to 0.1 with molecules oriented quite randomly after high temperature simulations on the system. There was no evidence of columnar or orientational ordering for this molecule.

Circum-coronene
The results for this molecule were quite different from those observed for the two smaller molecules. Plots of energy and order parameter are presented as a function of temperature in Figure 3. The energy was seen to increase with temperature; a very well defined phase transition could clearly be seen at T=0.2 that reached completion ~0.4. The system remained fairly stable at this energy and order parameter for T=2.0.
At P = 5, the structures obtained in simulations T = 0.10 to 0.18 were consistent with the ground state solid structure of the material, the only major difference being a few void spaces. From T = 0.2 to 0.28, an intermediate phase can be identified, followed by an isotropic liquid for T ≥ 0.30. Between phase at T = 0.18 and T = 0.2 there was little change in system energy; however, temperatures above this point showed a marked increase in system energy proportional to temperature. This change in slope suggests that a phase transformation may have occurred at this point. Further evidence of a phase change is the sudden drop in order parameter (from S eq = 0.995 to 0.967), as well as an increase in the order parameter's standard deviation by an order of magnitude (from σ S = 0.001 to 0.013). A discontinuity in an order parameter, with continuous system energy, is another feature typical of second-order phase transitions.
As indicated by the order parameter, the transition was associated with a decrease in the order of the molecular packing, as some molecules break away from the columns and into the voids between (see Figure 4). This behaviour represents the onset of fluidity within the system, suggesting that this phase represents a columnar-type mesophase. A clear columnar phase can be seen in Figure 4b.

Dodeca-benzo-circumcoronene
This regular molecule has the shape of a disc, quite large in size and appears to satisfy basic criteria for the formation of a mesophase. The evolution of order parameter as a function of scaled temperature is shown in Figure 5. Four distinct regions with changing order parameter can clearly be identified. Energy plots also showed a similar trend indicating various phase changes in the system. Figure 6 shows the centre of mass of molecule in the ground state and their initial relaxation as the temperature is increased slightly to 0.3 (7a). Although columns have tilted significantly, orientational order is being maintained to a great extent. As the temperature is increased further to T=1.1, columns are getting disturbed and molecular disorientation is setting in. The orientational order is still being maintained to a certain extent. At high temperature (T=2.0), system is completely disordered. Both centre of mass as well as molecular orientation is completely destroyed.

Representative pitch molecule
Two sets of studies were carried out on the representative pitch molecule. In one case, the aliphatic ring was treated as being distinct from the aromatic ring with additional hydrogen atoms being out of plane of the main molecule. Secondly, the aliphatic ring was treated as an aromatic ring resulting in a planar molecule. High temperature simulation results for both molecules are given in Figure 7. Three well-defined phases can be seen clearly for the nonplanar molecule while using NPT ensemble. These transitions were however not very well defined for constant volume simulations. This is to be expected as some of these transitions were accompanied by a volume change and could be suppressed to an extent in a fixed volume simulation.
For the non-planar molecule, three distinct discontinuities, suggesting possible phase transitions can be seen. The discontinuities in the energy plot also resulted in corresponding breaks in the orientational order parameter. The planar molecule showed only two  structural transformations for the range of temperatures and pressures investigated. Regime I in (b) appears to be quite similar to the regime I observed for the non-planar molecule. But the transition temperature (~0.6) in this case is much higher than the corresponding temperature (~ 0.27) for the non-planar molecule. Similarly the transition to regime II has also shifted to a much higher temperature (~0.95) as compared to the ~ 0.55 for the nonplanar molecule.
In Figure 8, we have plotted the positions of centres of mass, directors and atomic positions at T=1.15 for both molecules. These results point towards the important role played by the molecular structure wherein minor changes in structure can result in significant changes in stacking behaviour and structural transition points.

Conclusions
An extensive investigation has been reported on a range of molecules to determine key factors that could affect the formation of carbonaceous mesophase. The effect of molecular weight is a critical factor controlling mesophase formation from pitch precursors. However, in practical systems, molecular weights are not constant; polymerisation processes during heat treatment can cause molecular structures to change. Usually, in thermo-tropic liquid crystals, the mesophase forms upon heating of a solid or cooling of a liquid [31]: solid → mesophase → liquid. However, mesophase formation during the carbonisation of pitches usually occurs by a more complex route: solid pitch → isotropic liquid → mesophase → solid semi-coke. The three transitions in this process chain can be described as (a) Melting (or softening) -this effectively represents a reversible reduction in viscosity, (b) Mesophase formation -this step is characterised by an irreversible increase in viscosity and the formation of an anisotropic liquid (mesophase), (c) Carbonisation -this transition signifies the fusing of the material into a brittle coke The mesophase formation in carbon systems depends on more than simple thermo-tropic behaviour. At low temperature heat treatment, the pitch simply melts entering the isotropic liquid phase region. Once the temperature reaches a critical point, condensation and polymerisation begins, causing an increase in molecular weight which, in turn, makes the molecules more meso-genic. After sufficient polymerisation, the material enters the mesophase region. Further heat treatment increases the extent of polymerisation and cross-linking, fusing the material into a semi-coke.
Based on these experimental observations, it would be expected that the atomistic model as a representative of PAH systems, would show little or no mesophase formation for molecules of low molecular weight, representing the isotropic pitch melt at low temperatures (300-400 Da). Secondly, a strong tendency exists for mesophase formation for molecules of higher molecular weight, representing the partially polymerised pitch melt at higher temperatures (400-1000 Da) The absence of mesophase formation in coronene and hexabenzocoronene in atomistic simulations is in good agreement with these results. However, there was no evidence of nematic mesophase behaviour and only limited evidence of columnar mesophase behaviour in circumcoronene. These molecules have molecular weights of 300, 523 and 667 Da respectively. Since coronene lies within the molecular weight range for non-polymerised isotropic pitches, its inability to form a mesophase was to be expected. Other two molecules lie well within the higher weight range, where mesophase formation is to be expected.
Dodeca-benzo-circumcoronene, on the other hand, clearly showed nematic mesophase behaviour over a large temperature range. This molecule has a molecular weight of 1033 Da, placing it at the upper end of the stated range for meso-genic pitches. Thus, the model is in agreement with the experimental observations of pitches. Significantly different results were obtained for irregular, high-MW PAH (720 Daltons) representative of those in pitch. The simulations show columnar phases at low temperatures and the gradual loss of long-range columnar order with heating.
It is important to note that typical mesophase-forming pitches are complex mixtures; interaction between molecules of different sizes is a complex phenomenon. It is difficult to separate the effects of average molecular weight from the effects of molecular weight distribution and molecular shape. Thus, the study of the effects of molecular weight presented here should be combined with similar studies of the effects of molecular weight distribution and molecular shape; once this is achieved, it may be possible to create testable, quantitative predictions.