Phase Behaviors of Ionic Liquids Heating from Different Crystal Polymorphs toward the Same Smectic-A Ionic Liquid Crystal by Molecular Dynamics Simulation

Five distinct crystal structures, based on experimental data or constructed manually, of ionic liquid [C14Mim][NO3] were heated in NPT molecular dynamics simulations under the same pressure such that they melted into the liquid crystal (LC) phase and then into the liquid phase. It was found that the more entropy-favored structure had a higher solid-LC transition temperature: Before the transition into the LC, all systems had to go through a metastable state with the side chains almost perpendicular to the polar layers. All those crystals finally melted into the same smectic-A LC structure irrelevant of the initial crystal structure.


Introduction
Ionic liquid crystals (ILCs) are the intersections of ionic liquids (ILs) and liquid crystals (LCs), with the features of ILs being comprised of pure ions and of LCs having order(s) in one or two dimensions, whose application prospects [1,2] have been appreciated in the most recent twenty years.
Long-chain imidazolium-based ILs [C n Mim]X, with a proper combination of alkyl-chain length n and anion type X, are well known to be able to form LC structures and have been intensively studied. The combinations have been experimentally determined to be n ≥ 10 for (FH) 2 F, n ≥ 12 for NO 3 and BF 4 , n ≥ 14 for PF 6 , etc. [3]. Their thermotropic phase behavior, as an essential physical property, has been investigated by both experiments [3][4][5][6][7][8][9][10][11][12] and simulations [13][14][15][16][17]. In experiments, IL samples usually undergo several heating-cooling cycles, heated from a crystal phase to reach LC and/or liquid phases and then cooled down to recrystallize, during which crystal polymorphism is often reported, which is potentially of great importance in materials science and industry because physiochemical properties of polymorphs, such as thermal stability [6,9,18], ion conductivity [9,12], etc., can significantly affect productions and applications of ILs and ILCs. Due to the complexity of IL and ILC systems, however, the microscopic details of the phase behaviors related to IL polymorphs can only be speculative solely based on experimental data, in which molecular dynamics (MD) simulations can play a vital role [14,18,19].
From a simulation viewpoint, there are two convenient ways to obtain ILC structures. One is from heating a crystal structure [13][14][15][20][21][22], and another is through simulated annealing processes [17,23]. For the first way, the initial crystal structures are either handmade or from experimental data and may differ from one work to another, exhibiting the existence of crystal polymorphs. If the unit cells 1 Figure 1. Schematics of (a) enantiotropic and (b) monotropic systems. CrI and CrII are two crystal polymorphs. We propose that the ionic liquid (IL) polymorphism in this work belongs to the former, but the transition from CrI to CrII did not occur during our molecular dynamics (MD) simulations.

Simulation Setup
As we did in our previous works [14], we adopted an all-atom model based on an AMBER force field [34]. To obtain partial charges, we first performed an ab initio optimization at the mp2/6-31g* level with a Gaussian package [35], and then produced charges by fitting the electrostatic potential in space according to the Merz-Singh-Kollman scheme [36,37] at the same level. The partial charges were finally determined by the RESP utility [38] provided in the AMBER package. The MD simulations were run with the GROMACS 5.1.4 package [39]. The systems were coupled with an anisotropic or semi-isotropic Parrinello-Rahman barostat [40] to keep the shape of the simulation box flexible and properly commensurate with the crystal and LC structures. Each component of the reference pressure was set to 1 bar, and the compressibility was 4.5 × 10 −5 bar −1 . The temperature was controlled by a Nosé-Hoover thermostat [41], and a periodic boundary condition (PBC) was applied to all three dimensions. Electrostatic interactions were computed with the Particle-Mesh Ewald method [42] with a real-space cutoff of 12 Å, the same as the cutoff for van der Waals (VDW). VMD software [43] was used when visualizing simulation trajectories and when constructing and modifying crystal structures.
In the production simulations, temperatures were increased stepwise with a step of 25 K (50 K for the larger-size double bilayer (DB) system described below), and 5 K near the transition points. For the crystal phase, whose fluctuations were small, the systems were sampled within 10 ns. When the temperature approached the melting point and the systems became LCs or liquid structures whose fluctuations became large, the systems were sampled for more than 10 ns and up to 40 ns when necessary.

Preparation of Initial Configurations
The first handmade crystal structure, denoted by HM1 hereafter, was obtained by elongating the alkyl chains of the manually constructed [C 12 Mim][NO 3 ] structure we presented in a previous work [14]. The second handmade one (HM2) was newly constructed by greatly enlarging the tilt angle of the alkyl chains in HM1, so the layer spacing became much narrower. The other three structures were obtained by modifying the experimentally determined structures by replacing the anions with nitrates and removing the water molecules wherever they presented. DB represents the double bilayer structure modified from Reference [28], whereas EM1 and EM2 refer to the bilayer structures modified from References [28,29], respectively. The simulation boxes were made up by copying unit cells along the lattice vectors. The lattice vector a was always along the x axis, and b was in the x-y plane, so the polar bilayers were parallel to the x-y plane. The unit cells of HM1 and HM2 with 4 ion pairs were copied 6 × 6 × 3 times along lattice vectors a, b, and c, so the simulation boxes totally contained 432 ion pairs. A unit cell of a DB with 6 ion pairs was copied 4 × 4 × 4 times and 6 × 6 × 3 times to form two system sizes of 384 and 684 ion pairs, respectively. EM1, having 2 ion pairs per unit cell, was copied 6 × 4 × 8 times, and EM2, with 4 ion pairs, was copied 4 × 3 × 8 times, both resulting in a simulation box with 384 ion pairs.
With the anisotropic barostat applied, the crystal lattices were equilibrated at 5 K for 200 ps. Series of 2-ns NPT simulations then followed to increase the temperature stepwise to 300 K for further equilibration. The systems were annealed to 200 K, at which the production simulations began. The obtained crystal data of the polymorphs at 200 K are listed in Table 1, and the unit cells are drawn in Figure 2. Although the anions were replaced by nitrates, the DB, EM1, and EM2 retained their essential structural features with the space groups unchanged. The conformation of cation was also preserved: The cations were straight in HM1, HM2, and EM1; straight in one layer and bend (head group perpendicular to side chain) in neighboring layers in the DB; and "crank-handle"-like in EM2. with 4 ion pairs were copied 6 × 6 × 3 times along lattice vectors a, b, and c, so the simulation boxes totally contained 432 ion pairs. A unit cell of a DB with 6 ion pairs was copied 4 × 4 × 4 times and 6 × 6 × 3 times to form two system sizes of 384 and 684 ion pairs, respectively. EM1, having 2 ion pairs per unit cell, was copied 6 × 4 × 8 times, and EM2, with 4 ion pairs, was copied 4 × 3 × 8 times, both resulting in a simulation box with 384 ion pairs. With the anisotropic barostat applied, the crystal lattices were equilibrated at 5 K for 200 ps. Series of 2-ns NPT simulations then followed to increase the temperature stepwise to 300 K for further equilibration. The systems were annealed to 200 K, at which the production simulations began. The obtained crystal data of the polymorphs at 200 K are listed in Table 1, and the unit cells are drawn in Figure 2. Although the anions were replaced by nitrates, the DB, EM1, and EM2 retained their essential structural features with the space groups unchanged. The conformation of cation was also preserved: The cations were straight in HM1, HM2, and EM1; straight in one layer and bend (head group perpendicular to side chain) in neighboring layers in the DB; and "crank-handle"-like in EM2.

Order Parameters
Different phases have different spatial orders, which can be quantified by carefully defined order parameters. Distinct phases can be identified by a combination of several order parameters. To quantify the translational order of crystal structures, the translational order parameter (TOP) of a configuration is defined as where vector r denotes the position of the αth site in the jth unit cell, vector G is the sum of the three basis vectors of reciprocal lattice, | . . . | means taking the modulus, and n t and N UC are the number of a certain type of sites in the unit cell and the number of unit cells in the simulation box, respectively. When the systems are in the SmA LC phase, which has no crystal lattices but a one-dimensional translational periodicity, the TOP is modified accordingly as [44] where z is the z-component of the position, d is the layer spacing, and N t is the total number of cationic headgroups and anions in the simulation box. The basis vectors of the reciprocal lattice are determined from the real-space lattice. One set of the basis vectors of the real-space lattice is obtained by dividing the box vectors by the numbers of copies of unit cells along them. One of the basis vectors in the set is fixed in GROMACS to be along the x axis, and the other two vectors are further selected among the allowed primitive vectors (linear combinations of the other two basis vectors in the set) so that τ takes the maximum. Similarly, for LCs, the layer spacing d is chosen by maximizing τ LC .
With a directional vector defined for a cation connecting the center of masses (COMs) of its head and tail groups (the CH 3 group at the end of alkyl chain), the orientational order parameter (OOP) of a configuration is defined to quantify the orientational order of the system, which is where P 2 denotes the third term of the Legendre polynomials, θ is the angle between the directional vector of a cation and the average directional vector of the simulation box, and . . . means to take the average among the cations.
To quantify the conformational order of alkyl chains, all the dihedral angles defined by three consecutive carbon-carbon bonds on the same chain are computed, based upon which the ratio of the trans-conformation, whose angle is between 150 • and 180 • , is calculated and used as the conformational order parameter, denoted by R t .
These order parameters, together with the energy and density of the systems, are taken into account when judging whether the system is in equilibrium or at least in local equilibrium. The averages and variances of these quantities over nanosecond time scales are computed when their values with respect to time have no evident drift.

Results
When heating the systems, all the initial crystal polymorphs underwent phase transitions from crystal to LC and then from LC to liquid. Solid-solid phase transitions could happen below the solid-LC transition temperature, and a metastable state always occurred during the solid-LC transition. Note that the solid-solid phase transitions were not from one of the five initial structures to another, but rather the crystal structures initiated from the five simulated systems before melting were always different from each other. We present the snapshots of these phases during the simulations of HM1 in Figure 3. The snapshots of the metastable state, LCs, and NSL for the other four structures are basically the same. In the crystal phase (Figure 3a), all the polymorphs shared the same feature that cationic headgroups and anions arranged in an ordered pattern, forming polar bilayers with the nonpolar region composed of all-trans interdigital alkyl chains in between (except that EM2 had a gauche configuration near the headgroup). This structural feature held for the whole crystal range despite the fact that the lattice constants differed from polymorph to polymorph and changed during the solid-solid phase transitions. The DB structure was exceptional in that the number densities of cationic headgroups and anions in neighboring polar bilayers were different (see the unit cell in Figure 2): The layer with straight cations had a number density twice as much as the layer with headgroups bent to the alkyl chains, resulting in a doubled periodicity. In the metastable state ( Figure 3b) between crystal and LC, which had an intermediate order, the alkyl chains turned to be more perpendicular to the polar bilayers and no longer in the straight all-trans style, but they were still interdigital with each other. The LC structure ( Figure 3c) was even more disordered than the metastable structure, although the structural layering still existed, whose alkyl chains were no longer interdigitally aligned due to large thermal energy, and thus the layer spacing increased accordingly. As the temperature further went up, the polar layers dispersed, and the system eventually changed into the NSL phase ( Figure 3d). Crystals 2019, 9, x FOR PEER REVIEW 6 of 14 cationic headgroups and anions in neighboring polar bilayers were different (see the unit cell in Figure 2): The layer with straight cations had a number density twice as much as the layer with headgroups bent to the alkyl chains, resulting in a doubled periodicity. In the metastable state ( Figure  3b) between crystal and LC, which had an intermediate order, the alkyl chains turned to be more perpendicular to the polar bilayers and no longer in the straight all-trans style, but they were still interdigital with each other. The LC structure ( Figure 3c) was even more disordered than the metastable structure, although the structural layering still existed, whose alkyl chains were no longer interdigitally aligned due to large thermal energy, and thus the layer spacing increased accordingly.
As the temperature further went up, the polar layers dispersed, and the system eventually changed into the NSL phase ( Figure 3d). The transition temperatures and corresponding potential energy changes are listed in Table 2. The details of the structures appearing in the table could be characterized by the density, layer spacing, and order parameters, as is shown below. The phase transitions in Table 2 corresponded to the kinks on the caloric curves shown in Figure 4. The PV term was in fact negligible (~10−2 kJ/mol) compared to the large potential energy since the interactions in our model were strong and the applied pressure (1 bar) was relatively small, so the potential energy changes basically equaled the enthalpy changes. Note that in EM2, there were two solid-solid phase transitions, whose potential energy changes were both negative because the initial structure obtained by modifying the reported structure of [C14Mim][PF6] and re-equilibrating below 300 K was not the most stable structure for [C14Mim][NO3], and when the temperature was adequately high, the structure could transit into the more stable structure with a lower potential energy. Moreover, the LC-liquid transitions ca. 655 K were unidentifiable only by looking at the caloric curves in Figure 4, since the potential energy differences were too small, which agreed well with the experimental results of similar ILCs [4][5][6][7]10,12]. The transition temperatures and corresponding potential energy changes are listed in Table 2. The details of the structures appearing in the table could be characterized by the density, layer spacing, and order parameters, as is shown below. The phase transitions in Table 2 corresponded to the kinks on the caloric curves shown in Figure 4. The PV term was in fact negligible (~10−2 kJ/mol) compared to the large potential energy since the interactions in our model were strong and the applied pressure (1 bar) was relatively small, so the potential energy changes basically equaled the enthalpy changes. Note that in EM2, there were two solid-solid phase transitions, whose potential energy changes were both negative because the initial structure obtained by modifying the reported structure of [C14Mim][PF6] and re-equilibrating below 300 K was not the most stable structure for [C14Mim][NO3], and when the temperature was adequately high, the structure could transit into the more stable structure with a lower potential energy. Moreover, the LC-liquid transitions ca. 655 K were unidentifiable only by looking at the caloric curves in Figure 4, since the potential energy differences were too small, which agreed well with the experimental results of similar ILCs [4][5][6][7]10,12]. cationic headgroups and anions in neighboring polar bilayers were different (see the unit cell in Figure 2): The layer with straight cations had a number density twice as much as the layer with headgroups bent to the alkyl chains, resulting in a doubled periodicity. In the metastable state ( Figure  3b) between crystal and LC, which had an intermediate order, the alkyl chains turned to be more perpendicular to the polar bilayers and no longer in the straight all-trans style, but they were still interdigital with each other. The LC structure (Figure 3c) was even more disordered than the metastable structure, although the structural layering still existed, whose alkyl chains were no longer interdigitally aligned due to large thermal energy, and thus the layer spacing increased accordingly.
As the temperature further went up, the polar layers dispersed, and the system eventually changed into the NSL phase ( Figure 3d). The transition temperatures and corresponding potential energy changes are listed in Table 2. The details of the structures appearing in the table could be characterized by the density, layer spacing, and order parameters, as is shown below. The phase transitions in Table 2 corresponded to the kinks on the caloric curves shown in Figure 4. The PV term was in fact negligible (~10−2 kJ/mol) compared to the large potential energy since the interactions in our model were strong and the applied pressure (1 bar) was relatively small, so the potential energy changes basically equaled the enthalpy changes. Note that in EM2, there were two solid-solid phase transitions, whose potential energy changes were both negative because the initial structure obtained by modifying the reported structure of [C14Mim][PF6] and re-equilibrating below 300 K was not the most stable structure for [C14Mim][NO3], and when the temperature was adequately high, the structure could transit into the more stable structure with a lower potential energy. Moreover, the LC-liquid transitions ca. 655 K were unidentifiable only by looking at the caloric curves in Figure 4, since the potential energy differences were too small, which agreed well with the experimental results of similar ILCs [4][5][6][7]10,12].  The structural differences between polymorphs and the structural changes during phase transitions could be illustrated by density and layer spacing, shown in Figure 5. The layer spacing d in Equation (2) of a crystal is the height of a unit cell perpendicular to the polar bilayers, and that of an LC is the one-dimensional periodicity along the ordered direction. The density reflected the compactness of packing, and no experimental density data for [C 14 Mim][NO 3 ] were available for comparison. Figure 5 shows that the crystal polymorphs had similar densities, but different layer spaces. Sharp changes in density and layer spacing took place when the solid-solid and solid-LC phase transitions happened. In contrast, for the LC-liquid transition, the density changes were very small, and the fluctuations of the estimated LC layer spacing increased with temperature and became divergent at the clearing point.   The structural differences between polymorphs and the structural changes during phase transitions could be illustrated by density and layer spacing, shown in Figure 5. The layer spacing d in Equation (2) of a crystal is the height of a unit cell perpendicular to the polar bilayers, and that of an LC is the one-dimensional periodicity along the ordered direction. The density reflected the compactness of packing, and no experimental density data for [C14Mim][NO3] were available for comparison. Figure 5 shows that the crystal polymorphs had similar densities, but different layer spaces. Sharp changes in density and layer spacing took place when the solid-solid and solid-LC phase transitions happened. In contrast, for the LC-liquid transition, the density changes were very small, and the fluctuations of the estimated LC layer spacing increased with temperature and became divergent at the clearing point. The crystal phase was quantitatively characterized by computing the translational order parameters (TOPs). The crystalline TOP computed for the COMs of cations according to Equation (1) are shown in Figure 6a. For the HM1 and HM2 structures with 432 ion pairs, the TOPs kept above ca. The crystal phase was quantitatively characterized by computing the translational order parameters (TOPs). The crystalline TOP computed for the COMs of cations according to Equation (1) are shown in Figure 6a. For the HM1 and HM2 structures with 432 ion pairs, the TOPs kept above ca. 0.95 until they melted into LCs, justifying that they were strict in crystalline arrangement throughout the crystal phase. For the DB, EM1, and EM2 with 384 ion pairs, the TOPs were also beyond 0.7 before melting, although the values were relatively smaller, and the fluctuations were large. However, a series of additional simulations of the DB system with 648 ion pairs showed larger TOPs and smaller fluctuations, as can be seen in Figure 6b, indicating a noticeable finite-size effect. Nevertheless, the TOP values for the system with 384 ion pairs were already large enough to demonstrate that the COMs of cations were basically at the lattice positions before the melting points.
Similarly, the SmA LC was quantified by one-dimensional LC TOP, defined in Equation (2), calculated for the COMs of cationic headgroups and the COMs of anions. The LC TOP is also shown in Figure 6a. Compared to the large values of crystalline TOPs, the LC TOP was only about 0.5, indicating that the LC phase of the IL was rather disordered and merely had a rough layered structure. 0.95 until they melted into LCs, justifying that they were strict in crystalline arrangement throughout the crystal phase. For the DB, EM1, and EM2 with 384 ion pairs, the TOPs were also beyond 0.7 before melting, although the values were relatively smaller, and the fluctuations were large. However, a series of additional simulations of the DB system with 648 ion pairs showed larger TOPs and smaller fluctuations, as can be seen in Figure 6b, indicating a noticeable finite-size effect. Nevertheless, the TOP values for the system with 384 ion pairs were already large enough to demonstrate that the COMs of cations were basically at the lattice positions before the melting points. Similarly, the SmA LC was quantified by one-dimensional LC TOP, defined in Equation (2), calculated for the COMs of cationic headgroups and the COMs of anions. The LC TOP is also shown in Figure 6a. Compared to the large values of crystalline TOPs, the LC TOP was only about 0.5, indicating that the LC phase of the IL was rather disordered and merely had a rough layered structure. The orientational order parameters (OOPs) and the number ratios of the trans conformation of the alkyl chains, defined in Equations (2) and (3), respectively, are shown in Figure 7. As indicated by the OOP values shown in Figure 7a, which were close to 1 at low temperatures, the cations in the unit cell were almost along the same direction. The OOPs decreased sharply to about 0.5 when the crystals turned to LCs, suggesting that, unlike the common rod-like LCs, the soft alkyl chains of ILCs were not well aligned. When the temperature increased to around 655 K, the OOPs decreased suddenly to 0.1~0.2, corresponding to the occurrence of the LC-liquid transitions. The OOP values did not go to zero probably due to the finite-size effect. The orientational order parameters (OOPs) and the number ratios of the trans conformation of the alkyl chains, defined in Equations (2) and (3), respectively, are shown in Figure 7. As indicated by the OOP values shown in Figure 7a, which were close to 1 at low temperatures, the cations in the unit cell were almost along the same direction. The OOPs decreased sharply to about 0.5 when the crystals turned to LCs, suggesting that, unlike the common rod-like LCs, the soft alkyl chains of ILCs were not well aligned. When the temperature increased to around 655 K, the OOPs decreased suddenly to 0.1~0.2, corresponding to the occurrence of the LC-liquid transitions. The OOP values did not go to zero probably due to the finite-size effect. 0.95 until they melted into LCs, justifying that they were strict in crystalline arrangement throughout the crystal phase. For the DB, EM1, and EM2 with 384 ion pairs, the TOPs were also beyond 0.7 before melting, although the values were relatively smaller, and the fluctuations were large. However, a series of additional simulations of the DB system with 648 ion pairs showed larger TOPs and smaller fluctuations, as can be seen in Figure 6b, indicating a noticeable finite-size effect. Nevertheless, the TOP values for the system with 384 ion pairs were already large enough to demonstrate that the COMs of cations were basically at the lattice positions before the melting points. Similarly, the SmA LC was quantified by one-dimensional LC TOP, defined in Equation (2), calculated for the COMs of cationic headgroups and the COMs of anions. The LC TOP is also shown in Figure 6a. Compared to the large values of crystalline TOPs, the LC TOP was only about 0.5, indicating that the LC phase of the IL was rather disordered and merely had a rough layered structure. The orientational order parameters (OOPs) and the number ratios of the trans conformation of the alkyl chains, defined in Equations (2) and (3), respectively, are shown in Figure 7. As indicated by the OOP values shown in Figure 7a, which were close to 1 at low temperatures, the cations in the unit cell were almost along the same direction. The OOPs decreased sharply to about 0.5 when the crystals turned to LCs, suggesting that, unlike the common rod-like LCs, the soft alkyl chains of ILCs were not well aligned. When the temperature increased to around 655 K, the OOPs decreased suddenly to 0.1~0.2, corresponding to the occurrence of the LC-liquid transitions. The OOP values did not go to zero probably due to the finite-size effect. As for the number ratio of trans conformations, Figure 7b shows that the alkyl chains in the crystal phase, especially at low temperatures, were basically in the all-trans conformation. The ratio of EM2 crystal was an exception, and was about 0.909 (10/11), because of the gauche conformation composed of the first four carbons near the headgroup, i.e., the "crank-handle"-like conformation (see the unit cell shown in Figure 2). As the temperature increased, solid-solid phase transitions happened, and the alkyl chain conformations became less ordered, indicated by the decrease of the trans-conformation ratio. At the melting point, the decrease was more evident, because the breakdown of the interdigital structures allowed the alkyl chains to take more conformations without strong constraints between each other. On the contrary, the trans-conformation ratio remained unchanged at the clearing point, demonstrating that the conformations of alkyl chains in LC and in liquid were the same.

Discussion
It should be noted that, as a common problem for MD simulations, the transition temperatures were different from those obtained from experiments [3] because of hysteresis. It is a knotty problem to acquire transition temperatures in MD simulations comparable to experimental results for complex systems such as ILs, even with the aid of some free-energy-involved methods [45][46][47]. Nevertheless, the transition temperatures for all systems were systematically shifted to one side in the MD simulations, which did not hinder the qualitative comparison with each other.
There are two types of crystal polymorphisms, which are schematically shown in Figure 1. For the enantiotropic system, solid-solid phase transitions from one to another take place below the melting points, but not for the monotropic system. We suggest that the IL system in this work is an enantiotropic system. Note that the free energy at zero temperature equals the potential energy, and the slopes of the curves in Figure 1 equal the opposite of entropies. If CrI does not turn to CrII at the intersection of their free energy curves in Figure 1a in a finite-time simulation, it turns to LC at a lower temperature with a smaller entropy near the transition points, and CrII turns to LC at a higher temperature with a larger entropy, which coincides with the simulation data of the IL system here: Generally, the IL system with a lower density (Figure 5a) and smaller order parameters (Figure 7), indicating a more disordered structure with a larger entropy, tended to have a higher solid-LC transition point.
In our MD simulations, CrI and CrII could refer to any two of the five polymorphs, and we did not see the solid-solid phase transition from CrI to CrII, which should thermodynamically happen according to the diagram in Figure 1a. The reason was that the energy barrier between the two crystal structures was so high that this was almost impossible to occur during very limited MD simulation time and even experimental time [7]. Moreover, the type of anion may also have affected the occurrence of the solid-solid phase transition, as it has been reported to happen in ILs with BF 4 [12], but not with other anions. After all, the inhibition of transitions between different crystal structures leads to the counterintuitive phenomenon that more entropy-favored systems have higher solid-LC transition temperatures.
The processes of the solid-LC transitions were studied in detail. First, for the metastable states occurring during this transition, we could see by eye that the alkyl chains became upright. We therefore calculated the orientation of the cations to quantify this phenomenon and show the results in Figure 8.
As shown in Figure 8a, during the solid-LC transition, the system first turned to a metastable structure whose lifetime could be very short, and the alkyl chains became perpendicular to the polar layers, before it melted into the SmA LC, except the case of the DB, when the alkyl chains were already perpendicular in the crystal structure. Take the example of EM2, whose metastable state had a relatively long lifetime. From Figure 8b, in the crystal phase, both the polar and azimuth angles of alkyl chains were restricted in a narrow range. In the metastable state, the polar angles were still restricted to ca. 10 • , but the azimuth angles were evenly distributed in all directions, demonstrating that the alkyl chains were on average perpendicular to the polar layers. In the SmA LC, the polar angles increased to ca. 30 • with a broader distribution because the alkyl chains were no longer interdigital and bent. From a thermodynamic viewpoint, the distributions of alkyl-chain directions in metastable states were more symmetric than those in crystals, and the trans-conformation ratio of alkyl chains in metastable states was also smaller, leading to a larger entropy than crystal, but still smaller than LC. On the other hand, the alkyl chains in metastable state were still in the interdigital style, so the VDW potential between the alkyl chains, and thus the potential energy of the system, did not increase a lot. Together, the structure of metastable states favored entropy more without a large structural change costing the potential energy too much, so it was a reasonable intermediate structure along the solid-LC transition pathway. restricted to ca. 10°, but the azimuth angles were evenly distributed in all directions, demonstrating that the alkyl chains were on average perpendicular to the polar layers. In the SmA LC, the polar angles increased to ca. 30° with a broader distribution because the alkyl chains were no longer interdigital and bent. From a thermodynamic viewpoint, the distributions of alkyl-chain directions in metastable states were more symmetric than those in crystals, and the trans-conformation ratio of alkyl chains in metastable states was also smaller, leading to a larger entropy than crystal, but still smaller than LC. On the other hand, the alkyl chains in metastable state were still in the interdigital style, so the VDW potential between the alkyl chains, and thus the potential energy of the system, did not increase a lot. Together, the structure of metastable states favored entropy more without a large structural change costing the potential energy too much, so it was a reasonable intermediate structure along the solid-LC transition pathway. All the above quantities clearly show that, after the solid-LC transitions, all these crystal polymorphs transformed into the same SmA ILC. The curves for those quantities, which were separated from each other in the crystal ranges, overlapped very well after the solid-LC transition. This was true even for the DB structure with alternate anion densities in polar layers, as shown in Figure 9. When the crystal phase transitioned to the LC phase, the layer spacing increased, and the headgroups and anions moved across the nonpolar regions to rearrange, leading to nonzero densities of ions in nonpolar regions and then uniform polar layers with the same density. In the larger-size DB system (Figure 9b), this rearrangement completed shortly at the melting point of 530 K, with the anion densities in nonpolar regions coming back to zero, but it was more difficult for the smaller-size DB system (Figure 9a) to complete this process until 560 K. Consequently, the smaller DB system had all its quantity curves deviate from those for other systems in the ILC range from 530 K to 560 K, but it was a finite-size effect. All in all, we can conclude that heating a crystal IL can always produce the same ILC structure, and there is no need to consider the details of the initial crystal structure if we just want to have the ILC structure. It is probably the most thermodynamically stable structure in the temperature range from the melting point (up to 530 K) to the clearing point of about 655 K.
Based on our simulation results, we would like to remark on some experimental results in the literature. In the experiments of References [3][4][5][6][7][9][10][11][12] studying phase behaviors of ILs, the crystalline samples were measured during both heating and cooling. However, when we cooled the systems by MD simulation, they could only form SmA ILCs and dynamically go into a glassy state. Crystallization was practically almost impossible by MD simulation, let alone different crystal polymorphs, so we only presented the heating process. The comparisons between our simulation results and experimental results focused on layer spacing in the crystal phase or the ILC phase measured with X-ray scattering, as well as on the potential energy changes during phase transitions, All the above quantities clearly show that, after the solid-LC transitions, all these crystal polymorphs transformed into the same SmA ILC. The curves for those quantities, which were separated from each other in the crystal ranges, overlapped very well after the solid-LC transition. This was true even for the DB structure with alternate anion densities in polar layers, as shown in Figure 9. When the crystal phase transitioned to the LC phase, the layer spacing increased, and the headgroups and anions moved across the nonpolar regions to rearrange, leading to nonzero densities of ions in nonpolar regions and then uniform polar layers with the same density. In the larger-size DB system (Figure 9b), this rearrangement completed shortly at the melting point of 530 K, with the anion densities in nonpolar regions coming back to zero, but it was more difficult for the smaller-size DB system ( Figure 9a) to complete this process until 560 K. Consequently, the smaller DB system had all its quantity curves deviate from those for other systems in the ILC range from 530 K to 560 K, but it was a finite-size effect. All in all, we can conclude that heating a crystal IL can always produce the same ILC structure, and there is no need to consider the details of the initial crystal structure if we just want to have the ILC structure. It is probably the most thermodynamically stable structure in the temperature range from the melting point (up to 530 K) to the clearing point of about 655 K.
Based on our simulation results, we would like to remark on some experimental results in the literature. In the experiments of References [3][4][5][6][7][9][10][11][12] studying phase behaviors of ILs, the crystalline samples were measured during both heating and cooling. However, when we cooled the systems by MD simulation, they could only form SmA ILCs and dynamically go into a glassy state. Crystallization was practically almost impossible by MD simulation, let alone different crystal polymorphs, so we only presented the heating process. The comparisons between our simulation results and experimental results focused on layer spacing in the crystal phase or the ILC phase measured with X-ray scattering, as well as on the potential energy changes during phase transitions, more qualitatively, with the transition enthalpy measured with DSC. The comparison was limited to the species of [C 14 Mim]X with a size (ionic volume) of the anion X, similar to nitrate. The electronic structure of the anion and more detailed interactions of hydrogen bonding, anion-π, etc., which may affect the IL structure, especially the structure of polar layers [30], were not considered in this work.
To our knowledge, the solid-solid phase transitions of [C 14 Mim]X have been reported in three experimental works [6,9,12], in which the anions had similar sizes to nitrate. In Reference [6], solid-solid transitions were reported for [C 14 Mim] 2 [PdCl 4 ], which were irreversible for the low-temperature crystal grown from the solution and reversible for the high-temperature crystals occurring afterwards (along with the solid-LC transition). The layer spacing changed from 5.21 nm to 3.52 nm, and then to 3.72 nm for the crystal phases, and to 3.03 nm for the ILC phase. All these transitions had large transition enthalpies ranging from 16.1 to 18.8 kJ/mol, analogous to our observation of the DB structure melting into the LC. In Reference [9], solid-solid transitions were observed for [C 14 Mim][PF 6 ], where the layer spacing had a smaller increase and the transition enthalpies were reported to be much smaller than in Reference [6]. Although the EM2 structure was modified from the structure reported in this work, the transition observed in the experiment was more like the solid-solid transition of EM1 observed in our simulation by an increase in the degree of motion in the anion and alkyl chain after the transition indicated by order parameters. In Reference [12], the solid-solid transition changed the bilayer structure to a double or extended bilayer structure for [C 14 Mim] [BF 4 ]. The layer spacing of the bilayer structure was close to HM1, EM1, or EM2, and for the double or extended bilayer was close to the DB. Such a transition was not observed in our MD simulations.
Crystals 2019, 9, x FOR PEER REVIEW 11 of 14 more qualitatively, with the transition enthalpy measured with DSC. The comparison was limited to the species of [C14Mim]X with a size (ionic volume) of the anion X, similar to nitrate. The electronic structure of the anion and more detailed interactions of hydrogen bonding, anion-π, etc., which may affect the IL structure, especially the structure of polar layers [30], were not considered in this work. To our knowledge, the solid-solid phase transitions of [C14Mim]X have been reported in three experimental works [6,9,12], in which the anions had similar sizes to nitrate. In Reference [6], solidsolid transitions were reported for [C14Mim]2[PdCl4], which were irreversible for the low-temperature crystal grown from the solution and reversible for the high-temperature crystals occurring afterwards (along with the solid-LC transition). The layer spacing changed from 5.21 nm to 3.52 nm, and then to 3.72 nm for the crystal phases, and to 3.03 nm for the ILC phase. All these transitions had large transition enthalpies ranging from 16.1 to 18.8 kJ/mol, analogous to our observation of the DB structure melting into the LC. In Reference [9], solid-solid transitions were observed for [C14Mim][PF6], where the layer spacing had a smaller increase and the transition enthalpies were reported to be much smaller than in Reference [6]. Although the EM2 structure was modified from the structure reported in this work, the transition observed in the experiment was more like the solidsolid transition of EM1 observed in our simulation by an increase in the degree of motion in the anion and alkyl chain after the transition indicated by order parameters. In Reference [12], the solid-solid transition changed the bilayer structure to a double or extended bilayer structure for [C14Mim][BF4]. The layer spacing of the bilayer structure was close to HM1, EM1, or EM2, and for the double or extended bilayer was close to the DB. Such a transition was not observed in our MD simulations. Figure 9. Average number densities of nitrates along the z axis for (a) the DB structure with 384 ion pairs at 530 K in the crystal phase (black), at 530 K after the solid-LC transition (red), and at 560 K in the equilibrated LC phase (blue); and (b) for the DB structure with 648 ion pairs, in the crystal phase (black), after the solid-LC transition when the system was not well equilibrated (red), and in the equilibrated LC phase (blue), all at 530 K.
The metastable state was in accord with the transient state reported in Reference [11] and the low-temperature smectic phase in Reference [12], as described in our previous work [14]. When the metastable state melted into the SmA LC phase, the layer spacing of the ILC computed from simulations was close to the experimental results for the [C14Mim] salt with a similar-sized anion [6,7].
In experimental works [28], the DB structure has been proposed to explain the peak corresponding to spatial periods twice as large as those of bilayer structures in the direction perpendicular to the polar layers in X-ray scattering experiments, while another proposed structure has been the extended bilayer crystal structure [7,11,12], in which the alkyl chains are assumed to be in the end-to-end style, which was not observed in our MD simulations. The extended bilayer crystal structure seems unreasonable, since the area of the polar layers is incompatible with the cross-section of the nonpolar region if alkyl chains are in the end-to-end style, and thus the alkyl chains must bend Figure 9. Average number densities of nitrates along the z axis for (a) the DB structure with 384 ion pairs at 530 K in the crystal phase (black), at 530 K after the solid-LC transition (red), and at 560 K in the equilibrated LC phase (blue); and (b) for the DB structure with 648 ion pairs, in the crystal phase (black), after the solid-LC transition when the system was not well equilibrated (red), and in the equilibrated LC phase (blue), all at 530 K.
The metastable state was in accord with the transient state reported in Reference [11] and the low-temperature smectic phase in Reference [12], as described in our previous work [14]. When the metastable state melted into the SmA LC phase, the layer spacing of the ILC computed from simulations was close to the experimental results for the [C 14 Mim] salt with a similar-sized anion [6,7].
In experimental works [28], the DB structure has been proposed to explain the peak corresponding to spatial periods twice as large as those of bilayer structures in the direction perpendicular to the polar layers in X-ray scattering experiments, while another proposed structure has been the extended bilayer crystal structure [7,11,12], in which the alkyl chains are assumed to be in the end-to-end style, which was not observed in our MD simulations. The extended bilayer crystal structure seems unreasonable, since the area of the polar layers is incompatible with the cross-section of the nonpolar region if alkyl chains are in the end-to-end style, and thus the alkyl chains must bend and cannot result in a doubled periodicity. On the other hand, the experimental inference of heating a DB structure to ILC and then recrystallizing to a bilayer structure [7] seems to be contradictory to the fact that we have described that the DB structure is probably thermodynamically more stable than the bilayer structure at higher temperatures. A possible explanation is that the growth of the DB structure is easy in solvent [6,7], but very difficult in ILC, and hence it has never been reported. Another possibility concerns Ostwald's rule [48], according to which the least stable polymorph, here the bilayer structure, crystallizes first and becomes more stable as the temperature decreases, leading to the disappearance of DBs.

Conclusions
The thermotropic phase behaviors of IL [C 14 Mim][NO 3 ] during the heating of five distinct crystal polymorphs were investigated by MD simulation. The simulations started from the initial crystal structures at 200 K and ended at the clearing point of ILC, about 655 K, when the systems were in the liquid phase. The phases and phase transitions in between were recognized by measuring the potential energy, densities, layer spacing of crystals or LCs, and several order parameters. By analyzing the simulation data, we could reach the following three points. First, the IL is a special enantiotropic system characterizing an entropy-favored polymorph tending to have a higher melting point. Second, there always exists a metastable state between the crystal and the SmA LC, which has the feature that the alkyl chains become perpendicular to the polar layers and have more disordered conformations. Third, the structure of the SmA LC is irrelevant to the thermal history (i.e., different crystal polymorphs melt into the same SmA structure). The first point deserves further studies and may serve as an example for soft matter that entropy is of great importance in thermal stabilities. The second point presents a common state that all crystal polymorphs must go through during the melting of crystal into SmA LC, and may help uncover the melting process of similar ILC-forming ILs. The third point reveals that the SmA LC structure is independent of the crystal structure at a lower temperature, which is practically valuable for the study of ILC by MD simulation.
Author Contributions: W.C. and Y.W. conceptualized the research, W.C. performed the simulation and analyzed the data, and W.C. and Y.W. discussed the results and wrote the paper.