Computational Study on Homogeneous Melting of Benzene Phase I

Molecular-dynamics simulations are used for examining the microscopic details of the homogeneous melting of benzene phase I. The equilibrium melting temperatures of our model were initially determined using the direct-coexistence method. Homogeneous melting at a higher temperature is achieved by heating a defectand surfacefree crystal. The temperature-dependent potential energy and lattice parameters do not indicate a premelting phase even under superheated conditions. Further, statistical analyses using induction times computed from 200 melting trajectories were conducted, denoting that the homogeneous melting of benzene occurs stochastically, and that there is no intermediate transient state between the crystal and liquid phases. Additionally, the critical nucleus size is estimated using the seeding approach, along with the local bond order parameter. We found that the large diffusive motion arising from defect migration or neighbor-molecule swapping is of little importance during nucleation. Instead, the orientational disorder activated using the flipping motion of the benzene plane results in the melting nucleus.


Introduction
Melting is of fundamental importance in science because of interest in the manner in which crystalline solids spontaneously lose their energetically favored order [1,2].In nature, melting is usually initiated heterogeneously on surfaces, and a solid material immediately transforms into a liquid that is maintained at higher than the material's melting temperature [3,4].However, melting homogeneously occurs inside a crystal if surface melting is suppressed, which can be experimentally achieved using coating [5,6] or laser [7][8][9][10][11][12] techniques.Homogeneous melting can be defined as the phase transition from solid to liquid in the surface-and defectfree crystal, in which there is no trigger to initiate disordering except thermal fluctuations, and nucleation is equally likely to occur anywhere.
Benzene is a renowned hexagonal, small, and simple molecule.Regardless, its polymorphism [29][30][31] and local liquid structure [32][33][34] have not been completely explored.Although benzene's phase-transition dynamics have briefly been investigated [35][36][37], interesting phase behavior may be observed in the melting process.In 1982, Tohji et al. examined the lattice parameters of benzene phase I (stable phase under ambient pressure) via powder X-ray diffraction, and claimed the existence of a premelting stage that can be observed at 10 K below melting point [36].Furthermore, they predicted that the orientational disorder phase, plastic crystal, appears as a transient state during the transition process from the crystal phase to the liquid phase.The premelting effect was later considered to be disproven after an investigation over the complete temperature range of 4-280 K [37].However, the molecular mechanism of the melting dynamics of benzene, including the possibility of occurrence of a transient plastic phase, has not yet been clarified.
In the current study, we perform MD simulations for examining the molecular mechanism of the homogeneous melting of benzene phase I.We determine the equilibrium melting temperatures by implementing the direct-coexistence method and the homogeneous-melting temperature by heating up a surfacefree crystal.Further, we perform statistical analyses of induction times to investigate which of linear or stepwise transition is dominant in benzene melting.We used the seeding approach with the local bond order parameter for estimating critical nucleus size.Finally, we examine the key molecular motion that is required to form the critical melting nucleus.

Model
The benzene molecule was modelled with full atomistic detail using CHARMM22 [38,39].Van Der Waals energy is calculated with a standard 12-6 Lennard-Jones potential and the electrostatic energy with a Coulombic potential.It was demonstrated that this force field reproduces the crystal structure of benzene phase I [35,40].Other studies have also claimed that partial charges affect the orientation of benzene dimer [41] and the phase diagram for the other polycyclic aromatic hydrocarbons [42].Intermolecular interactions are truncated at 1.20 nm.The Lennard-Jones parameters of cross-interactions are obtained using the Lorentz-Berthelot combination rules: ij = ( ii • jj ) 1/2 and σ ij = (σ ii + σ jj )/2.Long-range Coulombic interactions are evaluated using the particle-mesh Ewald algorithm [43].

Molecular-Dynamics Simulation
MD simulations were conducted using the GROMACS 2018.3 package [44], where equations of motion are integrated with the leapfrog algorithm using a time step of 2 fs.Temperature T for production runs is controlled using the Nosé-Hoover thermostat [45,46] with damping constants of 1.0 ps, while the Berendsen algorithm [47] is used for equilibration.Pressure p is anisotropically controlled using the Berendsen algorithm [47] with damping constants of 2.0 ps, and pressurewasis set to 1.0 atm in all simulations.Energy minimization was conducted using the steepest-descent method.Periodic boundary conditions were applied in all three directions.
To determine the homogeneous melting temperature, we heated the system at 10 K intervals and performed an N pT-MD run for 10 ns at each T.This heating rate was chosen to observe nucleation after plenty of induction time, and to perform statistical analyses by making two hundred melting trajectories.

Crystalline Structure
The structure of the benzene phase I crystal unit cell was obtained from the Cambridge Structural Database (refcode: BENZEN, Deposition Number: 1108749) [48].The orthorhombic unit cell contains four benzene molecules, and the space group is Pbca.Further, lattice parameters are a = 7.460, b = 9.666, and c = 7.034 Å at 270 K [49].Each benzene molecules has 12 neighbors, and their centers of mass are arranged in a pseudo-face-centred cubic (pseudo-FCC) lattice.We replicated the unit cell in three dimensions to manufacture systems for analyses.In our simulations, lattice axes a, b, and c correspond to the x, y and z axes, respectively.

Melting Temperature
We initially evaluated melting temperature T m of benzene phase I in our model using the direct-coexistence method [50].We developed two identical boxes of benzene crystal containing 1344 molecules and melted one of them by performing a short NVT-MD simulation at 400 K. Box size was 5.4 × 5.8 × 5.6 nm.We joined these crystal and liquid boxes by erasing the overlapped molecules and making a 5.4 × 12.0 × 5.6 nm rectangular simulation box composed of 2539 molecules, named 'short system' (Figure 1A).We used the periodic boundary condition so that the liquid phase was in contact with the solid on both sides, although we do not show the imaginary cells.To relax the solid-liquid interface, we performed an NVT-MD run for 100 ps by freezing the solid phase, and the subsequent N pT-MD ran for 100 ps without freezing, where both simulations were conducted at 200 K. Further, we performed N pT-MD runs at 260, 265, and 270 K. Figure 1D depicts the time evolution of the average potential energy per molecule for these runs.The liquid has higher potential energy than the solid due to the disordered structure of the liquid, so we used the average potential energy to monitor which of melting or freezing occurs at a given temperature [50,51].At T m , the crystal and liquid phases exhibited the same stability in free energy, and their volumes did not change.On the other hand, if T > T m , the volume of liquid increased by the melting of the solid phase, so that the average potential energy increased.In contrast, if T < T m , potential energy decreased.Figure 1D shows that potential energy increased at 265 and 270 K, whereas it decreased at 260 K.These results indicate that 260 K is lower than T m , while 265 K is higher than T m , so that T m is somewhere between 260 and 265 K. Thus, we roughly assume that T m is 262.5 ± 2.5 K for the short system.
The presence of the interface might complicate the performance of N pT simulations [52].To assess the accuracy of our estimation for T m , we performed the same analyses for two other different cell sizes, medium (N = 5078) and long (N = 7617) systems.The initial configurations were prepared in the same manner as for the short system.The system size of medium and long systems was 5.4 × 23.0 × 5.6 nm, and 5.4 × 34.0 × 5.6 nm, respectively.These two systems had the same area of solid-liquid interface as the short system, and length in the direction perpendicular to the interface is two and three times longer, respectively (Figure 1B,C).In other words, these three systems have different size ratios of liquid-solid interface with respect to the whole simulation cell.The time evolution of potential energy for the medium and long systems indicates that T m is 265 ± 5 K (Figure 1E,F), which error range includes the T m estimated for the short system.Thus, we suppose that T m was 265 ± 5 K for our computational model, which is approximately 14 K lower than the experimental estimation of 278.7 K [53].
In a homogeneous regime, crystal structure can be superheated to a metastable state; its actual melting temperature is considerably higher than equilibrium melting temperature T m [5,10,15,20,50].To determine homogeneous melting temperature, we made a crystal structure of N = 7280, with a box size of 9.6 × 9.4 × 9.5 nm, and equilibrated the system at 240 K.Then, we heated the system in increments of 10 K, performing an N pT-MD simulation for 10 ns at each T. We observed that melting occurs at 330 K as soon as T is increased from 320 K. Figure 2A depicts that the average potential energy per molecule jumps between T = 320 and 330 K owing to the first-order phase transition, whereas it continuously increases at other temperatures.This result indicates that the limit of superheating achieved by our heating rate (1.0 K/ns) is 330 K, which is approximately 70 K higher than T m .Furthermore, we confirmed that the lattice parameters of the crystal (a, b, and c) change continuously at lower than 320 K (Figure 2B).Craven et al. demonstrated the smooth variation of the lattice parameters below T m using neutron-powder diffraction [37].In our result, smooth variation was observed until 320 K. Thus, we could conclude there was no hint of premelting, even for a superheated condition.

Induction Time
We investigated the characteristics of melting dynamics in terms of induction (waiting) time.Below the superheating limit, there is a free-energy barrier between the crystal and liquid phases, and the system was observed to remain at a superheated crystal phase until a melting pathway can be obtained.We obtained 100 independent melting trajectories at 323 and 324 K by providing a different set of velocities to the initial configuration, equilibrated at 320 K. Figure 3A displays the time evolution of the potential energy per molecule for three typical melting trajectories obtained at 323 K.During the induction period, the potential energy fluctuates; however, its mean value does not change.The blue curve shows a small jump of potential energy before melting starts, which arises from a longevity melting nucleus.However, this behavior was not observed in the remaining 199 melting trajectories.Thus, we did not give special weight to the jump.After several hundred ps, the ordered crystal structure spontaneously collapses with increasing potential energy.If there is a transient metastable state between the crystal and liquid phases, potential energy presents a stepwise change owing to the Ostwald step rule [54].However, the homogeneous melting of benzene in the 200 melting trajectories observed in our study was always completed without any steps.To further confirm that the homogeneous melting of benzene is governed by a stochastic transition between two states, we investigate the statistics of induction time τ.If the phase transition is random between the initial and final states, and exhibits a short memory time with a constant probability of occurrence per unit time (τ −1 0 ), the probability distribution of τ can be expressed by τ −1 0 exp(−τ/τ 0 ); further, average τ is τ = τ 0 [55].In addition, nucleation rate J is estimated by J = τ −1 0 V −1 , where V is the volume of the simulation box in the crystal state.However, if the transition is complicated, such as a multistep transition, the probability distribution of τ obeys a different function [56], and the computed values of τ 0 and τ do not agree.
We define induction time τ as the time period from the beginning of the MD run at 323 K (or 324 K) to the instant at which the potential energy is midway between the mean values of the crystal and liquid phases.For example, the τ of the melting trajectories in Figure 3A is 0.9, 5.5, and 8.3 ns. Figure 3B depicts the probability distribution of τs at 323 and 324 K.We observed that the probability distribution can be fitted well using the aforementioned exponential function (τ −1 0 exp(−τ/τ 0 )) at both temperatures even though distribution is strongly dependent on temperature.The estimated values of τ 0 s agree well with the values of τ , which are directly computed from MD runs; it was observed that τ 0 = 3.2 ns and τ = 3.4 ns at 323 K, and that τ 0 = 1.1 ns and τ = 1.2 ns at 324 K.We also estimated the nucleation rate to be J = 3.2 × 10 26 cm −3 s −1 at 323 K and 9.3 × 10 26 cm −3 s −1 at 324 K.These results indicate that the homogeneous melting of benzene can be described by a classical two-state transition model [57], in which the metastable crystal phase stochastically transforms to the stable liquid phase.

Nucleation
Here, we investigate the spontaneous formation of liquid nucleus inside the crystal.To distinguish the disordered structure from the structure of the benzene crystal with a pseudo-FCC lattice, we computed averaged bond order parameter q 6 for each benzene molecule using their centers of mass [58].q 6 captures local symmetry.Two benzene molecules are considered to be adjacent if the distance between their centers of mass is shorter than 0.69 nm, where the radial-distribution function of the crystal structures at 320 K exhibits a second minimum; further, the number of neighbors is 12 molecules.In Figure 4, we compare the probability distribution of q 6 for structures of bulk crystal at 320 K and bulk liquid at 330 K. Benzene-crystal distribution has a peak at q6 = 0.49, which is the same as that for a perfect FCC crystal [58].Further, there is no overlap between distributions.Therefore, q 6 is expected to identify the disordered benzene structure that is formed in the crystal.Here, we define a benzene molecule as liquid if its q 6 is less than 0.4.  .Probability distributions ofthe average bond order parameter q 6 for a bulk crystal at 320 K and a bulk liquid at 330 K.
We focus on one of the typical melting trajectories at 323 K, where melting occurs after sufficient induction time.The analysed trajectory is the same as that of the green line in Figure 3A. Figure 5A presents the time evolution of the total number of liquid benzene molecules n t and the largest cluster size n c , where two benzene molecules are considered to form the same cluster if they are adjacent.In Figure 5A, the spikes in n t and n c demonstrate that the melting nucleus, which can be considered to be a defect if the size is small, repeatedly appears and disappears during induction time.At this temperature, liquid benzene always exists in the crystal, and average n t and n c before melting (<8.2 ns) were 16.1 and 8.4 molecules in our system (N = 7280), respectively.When n t is small, there are clusters so that n c and n t are different (Figure 6A-C).When melting begins, one of them grows into the melting nucleus, and the remaining small clusters are subsequently merged into the largest (Figure 6D). Figure 6A-C shows that their birthplaces are not fixed at a specific position but randomly scattered in space.At approximately 8.25 ns, n c is equal to more than 200 molecules; further, the melting nucleus grows rapidly and transforms to bulk liquid within 150 ps (the inset of Figure 5A).These results imply that critical nucleus size is likely to be smaller than 200 molecules, beyond which the melting process is shifted from nucleation to nucleus growth.

Critical Nucleus Size
To determine the critical nucleus size, we adopted the seeding approach [59,60].A seed for a melting nucleus could be obtained from our melting trajectory.Further, 10 independent, 20 ps-long N pT-MD runs are initiated from the seed owing to the generation of different velocities.Figure 7A depicts the time evolution of n c for a swarm of trajectories beginning from the melting nucleus of n c = 172, corresponding to Structure C in Figure 6.The results denote that the melting nucleus grows in some trajectories, whereas it shrinks in other trajectories.The average value, n c , did not significantly change (see blue line in Figure 7A), thereby indicating that the system is close to the hill-top on the free-energy landscape.In comparison, the melting nucleus grows in the majority of the trajectories when the trajectory begins from the melting nucleus of n c = 210 (Figure 7B).Here, the initial nucleus size is already greater than the critical nucleus size.Further, we quantified average growth speed d n c /dt for an ensemble of seed nuclei ranging from n c = 60 to 220, randomly sampled from the obtained melting MD trajectories.Figure 7C depicts that d n c /dt is directly proportional to the initial n c even though its large dispersion indicates that n c is not the only order parameter defining the melting progress.The linear fit of the dataset shows critical nucleus size n * c with d n c /dt = 0 being approximately 110 molecules.This size is consistent with the observation that the fluctuation of n c rarely results in more than 110 molecules during induction time (Figure 5A).

Diffusive Motion
What is the key motion that is required for producing the critical melting nucleus?Studies on ice [20], metals [17,[21][22][23], Lennard-Jones fluids [24][25][26][27] and colloids [9] denoted that the precursor for homogeneous melting is the diffusion of molecules, which is activated by defects or co-operative motion.These precursors promote the swapping of molecules and change their lattice sites while the whole crystal structure remains intact.Further, we compute the time evolution of mean square displacement (MSD) with respect to the initial configuration, |r i (t) − r i (0)| 2 , where r i (t) is the coordinate vector of the centre of mass for the ith molecule at time t; further, bracket • • • denotes the average over all molecules.Figure 8A depicts that MSD does not increase during induction time, which is in contrast to that observed in a Lennard-Jones crystal [27].Furthermore, the projection against n c demonstrates that MSD increases after n c becomes larger than critical nucleus size n c * (Figure 8B).These results indicate that large diffusive motions induced by defect migration or neighbor-molecule swapping do not contribute to the formation of the melting nucleus.Furthermore, we compare two particular moments, III (t = 8.243 ns) and V (t = 6.740 ns), which are indicated by black arrows in Figure 8B.Although n c values at both moments were larger than n * c , the melting nucleus did not grow after Moment V. Further, we can observe that the MSD at Moment V remained small, indicating that the distance from the perfect crystal on the conformational space was not sufficiently large.Indeed, the system repairs the disordered structure and the crystal structure is restored by thermal fluctuation (see Figure 5).The difference between Moments III and V indicates that the quality of the melting nucleus, in addition to n c , should be taken into consideration for more accurately defining the melting progress.

Rotational Motion
We investigated the rotational motion of benzene crystal molecules.To characterize their orientation, we focused on two orthogonal unit vectors, vector q perpendicular to the molecular plane, and in-plane vector u that passes through two carbon atoms that are located on opposite sides of the molecule (Figure 9A).Vectors q and u denote the six-and twofold axes of the molecule, respectively.Figure 9B depicts self-reorientational correlation function (RCF) u i (t) • u i (0) for vector u at 240, 280, and 320 K.All functions decay in the order of a few tens of picoseconds, and the rotation of u becomes faster with increasing T. We evaluated relaxation time τ u by a fit of exp(−t/τ u ); values of τ u were 64, 22, and 9 ps at 240, 280, and 320 K, respectively.Even at the lowest temperature of 240 K (<T m ), relaxation time τ u = 64 ps was considerably shorter than heating time (10 ns long MD run at each temperature).Further, we analyzed the rotational behavior of one molecule.In Figure 9C,D, time evolution of angle θ u = acos{u i (t) • u i (0)} is plotted for a typical molecule at 240 and 320 K, respectively.The results depict that rotation occurs by discrete jumps between θ u = 0, 60, 120, and 180 degrees, and the frequency of the jumps increases with T owing to larger thermal energy.We later depict that the relaxation of q is hundreds of times slower than that of u, so that the reorientation of u is considered independent from that of q, occurring around the sixfold axis (vector q).In such a case, the conformation of the benzene molecule is identical between θ u = 0, 60, 120, and 180 degrees and does not distort the crystal structure.Because benzene molecules in a crystal frequently rotate around the sixfold axis, even at almost one hundred degrees lower than the homogeneous melting temperature, this motion is unlikely to influence the formation of the melting nucleus.RCF analyses depict that the relaxation of q was significantly slower than that of u (Figure 9E).At 240 and 280 K, the RCF does not decrease over a period of 10 ns, indicating the reorientation of q is negligible.However, the RCF decays in the order of a few ns at 320 K, which is close to the homogeneous melting temperature.We monitored the time evolution of angle θ q = acos{q i (t) • q i (0)} for each benzene molecule.The angle changed only a few times in 10 ns at 320 K. Figure 9F denotes one of the typical examples.The change between θ q = 0 and 180 degrees arises from the flipping motion of the molecule plane.If benzene molecules are completely flipped, the motion does not affect the crystal structure.However, we observed that benzene molecules sometimes fail to completely flip, thereby inducing orientational disorder.The benzene molecule is not spherical; the unfavored orientation dislocates the surrounding molecules, and the resulting asymmetric structure is detected by order parameter q6 .The overlap of the small melting nucleus, corresponding to Moment I in Figure 6, on the crystal structure is depicted in Figure 10.The snapshot clearly illustrates that the melting nucleus comprises orientational disorder, while the translational order of the centers of mass is maintained considerably well (we can observe the alignment along the vertical and horizontal directions).To further confirm the influence of the flipping motion on nucleation, we computed the average change of vector q at time t with respect to ∆t=10 ps before; ∆θ q (t) = ∑ N i acos{q i (t) • q i (t − ∆t)}/N.Then, we found that the correlation coefficient between ∆θ q (t) and the total number of liquid benzene n t was 0.51 at the induction time, indicating moderate correlation.Along with the analyses for diffusional and rotational modes, we can conclude that the flipping motion that was activated near the melting temperature triggered the formation of the melting nucleus.

Summary
In this study, the molecular mechanism of the homogeneous melting of benzene phase I was investigated using MD simulations.Further, we intended to understand the manner in which a nonspherical apolar organic crystal with a well-packed structure generated internal disorder.We demonstrated that melting occurs stochastically in the benzene crystal, and that there was no transient state between them at the superheated temperature at which our simulations have been performed.While defect migration and neighboring-molecule swapping were not observed before melting began, benzene molecules began to explicitly diffuse when the melting nucleus became larger than the critical size.In comparison, we observed that the flipping rotational motion is important for the formation of critical melting nucleus because the nucleus exhibits orientational disorder.In addition to elucidating the homogeneous melting of benzene phase I, our study may serve as the basis for understanding the phase-transition dynamics of other polycyclic aromatic hydrocarbons (such as perylene [61,62]) that possess a complicated free-energy landscape, including metastable crystals.

Figure 1 .Figure 2 .
Figure 1.Estimation of T m of benzene phase I in our computational model using the direct-coexistence method, performed for three different systems, (A) short, (B) medium, and (C) long.Initial configurations including liquid (right half) and solid (left half) phases, and time evolution of average potential energy (PE) per molecule in N pT-MD simulations at three temperatures (260, 265, and 270 K) for (D) short, (E) medium, and (F) long systems.

Figure 3 .
Figure 3. Induction time.(A) Time evolution of potential energy (PE) per molecule for the three melting trajectories obtained at 323 K; (B) probability distributions of induction times (τs) for melting at 323 and 324 K, where solid and dashed curves denote exponential functions τ −1 0 exp(−τ/τ 0 ) fitted to the distributions.

Figure 4
Figure 4. Probability distributions ofthe average bond order parameter q 6 for a bulk crystal at 320 K and a bulk liquid at 330 K.

Figure 5 .Figure 6 .
Figure 5.Time evolution of total number of liquid benzene molecules n t and largest liquid cluster size n c for (A) the whole induction time and (B) the short period just before melting occurs.Inset in (A) shows the entire melting process.Dashed line represents the size of critical nucleus n * c estimated by the seeding method (110 molecules).(A)(B)

Figure 7 .
Figure 7. Time evolution of nuclei in 10 independent trajectories, beginning from (A) n c = 172 and (B) 210.Initial seed in A corresponds to Snapshot C in Figure 6.Blue lines indicate average value n c .Red dashed lines are linear fits of n c , from which average growth speed d n c /dt can be calculated.(C) Plot of d n c /dt as a function of the initial n c of the seeds.Further, black dashed line indicates linear fit.

Figure 8 .
Figure 8. Mean square displacement (MSD), as computed for the same trajectory depicted in Figure 5, is plotted against (A) time and (B) largest cluster size n c .In (B), the dashed line indicates critical nucleus size n * c , and black arrows indicate Moments III and V that correspond to those observed in Figure 5

Figure 9 .
Figure 9. Analyses of rotational motions.(A) Illustration of two vectors u and q on a hexagonal benzene molecule.Reorientational correlation functions (RCFs) for vectors (B) u and (E) q at three temperatures.Time evolution of angles θ u in the unit of degree at (C) 240 K and (D) 320 K, and (F) angle θ q at 320 K, as computed for a typical benzene molecule.

Figure 10 .
Figure 10.Projection of the largest cluster (thick red) observed at Time I on the crystal structure (thin black).