The role of the dipole moment orientations in the crystallization tendency of the van der Waals liquids – molecular dynamics simulations

Computer simulations of model systems play a remarkable role in the contemporary studies of structural, dynamic and thermodynamic properties of supercooled liquids. However, the commonly employed model systems, i.e., simple-liquids, do not reflect the internal features of the real molecules, e.g., structural anisotropy and spatial distribution of charges, which might be crucial for the behavior of real materials. In this paper, we use the new model molecules of simple but anisotropic structure, to studies the effect of dipole moment orientation on the crystallization tendency. Our results indicate that proper orientation of the dipole moment could totally change the stability behavior of the system. Consequently, the exchange of a single atom within the molecule causing the change of dipole moment orientation might be crucial for controlling the crystallization tendency. Moreover, employing the classical nucleation theory, we explain the reason for this behavior.

the non-bonded interaction between atoms of different molecules, have been defined by OPLS all-atom force field 44 parameters provided for carbon atoms of the benzene ring. To keep simplicity of constructed model molecules as much as possible, we decided not to add hydrogen atoms. As a consequence, we obtained systems of rhombus-like molecules (RLM), which to some extent mirror the flat and asymmetric shape of the real molecules. Subsequently, we redefine the charges for appropriate atoms of RLM to ±0.125e, ±0.25e, or ±0.5e (e is an elementary charge) to create weak and strong μ oriented alongside the long and the short molecular axis. Finally, we obtained 5 different systems; the schemes of their constructions are presented in the Fig. 1. RLM of type A does not have a dipole moment, the RLM of type B possess the dipole moments oriented alongside the long molecular axis, whereas RLM of type C are characterized by the dipole moment of orientation along the short molecular axis. Moreover, molecules B1 and C1 are characterized by the same value of μ, which is two times smaller than μ for molecules B2 and C2.
Our investigation we began from the construction of perfect FCC lattice crystal, in which the 2048 RLM were inserted in the place of atoms. The standard simulations of the molecular dynamics have been performed by use of the GROMACS software 45,46 at conditions of constant temperature and the pressure controlled by the Nose-Hoover thermostat [47][48][49] and Martyna-Tuckerman-Tobias-Klein barostat 50,51 (p = 100 MPa) respectively. Each simulation run lasts for a relatively long time, i.e, 1 bilion of time-steps (dt = 0.001 ps). The first half of simulation was devoted for equilibration of the system, whereas the volumetric and dynamic data have been collected for the last 500 millions of time-steps. Then, we heated (ΔT = 5 K) the systems at constant pressure from the starting temperature equals 10 K. At this point we would like to mention that the initial FCC structure were proposed by use due to computational reasons. It protects system to undergo significant changes of the volume during the first steps of simulations at 10 K. Additionally, as one can see in Fig. 1a the initial structure forces systems to maintain some crystal structure, which immediately results that during heating process we observe the step increase in the volume for each system, which indicates on the melting of the crystal structure, see Fig. 1a where the temperature dependence of the molecular volume V mol = V/N (N denotes number of molecules within the system of volume, V) is presented. The heating was stopped at the temperature at least 35 K higher than the temperature of the detected drop in the density. During this stage one can observed that after equilibration at starting temperature the systems without μ, and with μ oriented Figure 1. The temperature dependences of molecular volumes during heating (a) and cooling (b) for all studied systems are presented. The star denotes the initial FCC configuration, which was subsequently equilibrated at T = 10 K and p = 100 MPa. In panel (b) the temperature at which the crystallization takes place is indicated. alongside the short molecular axis, i.e., systems A, C1, and C2 tend to form the more preferable (and denser) crystal structure than those resulted from initial FCC structure, whereas system B1 and B2 does not exhibit any phase transition to the another crystalline form. At this point we would like to put readers attention on the fact that some structural changes may occur also during equilibration process at 10 K. This fact enables us to expect that close to temperature of observed step increase in the volume analyzed systems exhibit thermodynamically stable crystal structures. Of course, those structure probably are not perfect. However, simplicity of designed molecules let us suspect that unit cell of their crystal form would not be complex and therefore the number of defects in comparison to the number of crystal cells is tiny. Consequently, we expect that influence of defects is negligible and therefore for our further theoretical analysis we will base on crystal structures registered during heating. Moreover similarly to the previous results for simple-liquids [19][20][21] , increase in the dipole moment, which implies gain in the attraction between molecules, causes that crystal phase become more stable, i.e., the temperature at which melting of the crystal structure is observed, become higher. Therefore, on the base of presented studies, we can confirm that this effect is independent on the orientation of the dipole moment.
After heating of the examined systems we started the main stage of our experiment, which is the isobaric cooling up to the starting temperature (ΔT = 5 K), see Fig. 1b. Interestingly, during cooling, the only one system exhibits the step decrease in the volume, which suggests the occurrence of the crystallization event. It is observed for system (C2), which is characterized by strong μ oriented alongside the short molecular axis. At this point, we would like stress that in order to verify whether this finding is repeatable, we performed three times cooling of systems C2 and B2 (for which temperatures of observed melting is the most nearing to temperature of observed melting of system C2) and all experiments ended identically.
Considering the results of systems cooling, which are comprised of various RLM we would like to put readers attention on certain key observation. Since the system C2 possesses higher value of μ than system C1, it is suspected that system C2 is characterize by the stronger attraction between molecules, and then, according to results of simple-liquid simulations, its tendency to crystallization should be smaller [19][20][21] . However, system C2 exhibits the sudden drop in volume indicating on crystallization whereas such behavior is not observed for system C1. Therefore, the straightforward relation between the gain in strength in the intermolecular attractions, which is induced by the value of the dipole moment, and tendency of the system to the formation of the crystalline structure is not evident for the investigated herein systems. Hence, behavior of RLM is in contrast to the recent results for the simple-liquids. As a consequence, the validity of discussed straightforward relationship cannot be expected for the real materials in general. Furthermore, this outcome suggests that the simulations of model systems, which are characterized by more complex architecture than the simple-liquids, e.g., RLM, could be a method for systematic and detailed studies of given intramolecular feature effect on structural, dynamic, and thermodynamic properties of the system. In this context it is worth to noting that the simplicity of RLM gives promising hope for theoretical solutions, which are not accessible for real materials (in this context we encourage readers to acquaint with our very recent paper where RLM are used to solving the important problem of the density scaling 43 ).
It is also worth noting that during cooling of systems A, B1, B2, and C1 the change in the slope of the temperature dependences of the volumes, which is a manifestation of the glass transition, is observed. It implies that RLM seem to be interesting candidates to study the glass transition as well, especially because the common one-component model systems usually easily crystallize.
Additionally, it should be noted that system C2 likely exhibits the highest trend to crystallization independently on the pressure conditions because, especially for simple molecules, increase in the pressure facilitates the formation of the crystal form within the liquid. Thus, also at lower pressures exclusively system C2 is suspected to occur liquid-crystal phase transition. However, this key observation of our studies as well as aforementioned advantages of RLM could be recognized as valid only if the step decrease in the volume observed for system C2 in Fig. 1b is indeed a manifestation of the discussed process. Given this remark, we analyzed the structure of the system C2 at crucial thermodynamic conditions. Already at the first glance a comparison between the simulation box snapshots at T = 135 K (, which is just before the expected crystallization event, T = 130 K), T = 125 K, and at the final temperature (T = 10 K), reveals evident differences, see insets of Fig. 2a-c. RLM of system C2 at T = 125 K and at the final temperature are visibly arranged in manner characteristic for the crystal structure, whereas any similar arrangement is not observed at T = 135 K. At this thermodynamic condition, the structure of system C2 seems to be totally disordered (like in the case of liquid), which indicates on the process of rapid structure ordering. Above visual conclusion is in accord with the observed changes in the shape of the radial distribution function (RDF) calculated at the mentioned thermodynamic conditions. It is clearly seen in Fig. 2b,c that the shapes of RDF exhibit the distinctive for the crystal structure the sharp peaks. On the other hand, the molecular configuration at T = 135 K is characterized by RDF of much softer shape, which is typical for liquid structure. Moreover, the values of RDF at this temperature are about two times smaller than at the final temperature. Thus, the changes in RDF also confirm that observed in Fig. 1b steep decrease in the volume for system C2 is indeed the manifestation of the crystallization process.
Another interesting point, which is worth to study, is the examination of the temperature behaviors of systems comprised of various RLM. However, the comparison between different shapes of RDF function is not the convenient and precise way for this analysis. Fortunately, it can be done using another standard tool for a characterization of the local/global structures which are the bond-orientational order parameters 52 . It is worth to note that these parameters are commonly used for identification of different crystalline phases [53][54][55][56] as well as to study of the melting transitions [57][58][59] . For this purpose, i.e., for the studies of supercooled liquids and glasses, q 6 and Q 6 are considered as the most prominent parameters, the idea of which could be briefly described as follows. The local structure around particle i could be quantified by a set of numbers where l = 6 for q 6 , r iĵ is a unit vector in the direction between particle i and one of its neighbours, j, which total number equals N(i), and , which is finally recognized as a global order parameter (Q 6 for = l 6). The temperature dependence of the global bond-orientational order parameter, Q 6 , between the centers of the molecules is presented in the inset of Fig. 3 for the system C2 and B2. For this study, we defined the nearest neighbor as the particle which is not farther than 1/10 of the simulation box from the considered particle. One can see that during cooling Q 6 increases in a consistent manner up to T = 135 K, i.e., up to the temperature below which the system C2 crystalizes. Then, Q 6 rapidly rises, which means immediate ordering of the structure -crystallization. After this process the structure of the system remains practically unchanged. Interestingly any discontinuity (jump) in the value of Q 6 is not detected for system B2. In this case Q 6 increases up to T = 80 K and then remains constant. This behavior indicates freezing of the liquid structure. It is worth to mention that RDF for system B2 at the final temperature still possesses characteristic for liquid structure shape (results not presented). Interestingly the freezing of the structure takes place at the temperature for which characteristic for the glass transition change of the slope of the temperature dependence on volume, is observed.
Since from the five examined systems only one crystallizes, it is interesting to check whether the occurrence of crystallization process could be explained within the classical nucleation theory (CNT), which is the most widespread theoretical description of the discussed process. According to the CNT the crystallization is preceded by the formulation of the nuclei of critical size, which subsequently grows. Therefore, the nucleation rate of critical nucleuses, J, is the crucial quantity determining the occurrence of the crystallization event. In the term of CNT, the nucleation rate is expressed according to the following formula 60 ,  In the insets, the snapshots of the simulation box at respective thermodynamic conditions are shown. Additionally, the scheme of RLM of system C2 is presented in Fig. 2a. www.nature.com/scientificreports www.nature.com/scientificreports/ zation). Before the determination of the Δg, it is usefull to consider that for all examined systems we set identical isobaric condition, and therefore used pressure could define as the reference pressure. Then the driving force become a function of only temperature, and consequently the equation for driving force proposed by Gutzow 61 with the redefined integration pathway 62 takes the following form g T S T dT where ΔS is the difference in the entropy between the liquid and the solid phases. ΔS was estimated using the obtained directly from the simulation-runs temperature dependence of the enthalpy, H T ( ) and the classical relation between the enthalpy and the entropy, T H S ( / ) p = ∂ ∂ , (considering the melting boundary condition, S H T / m m m Δ = Δ ). At this point it has to be stressed that calculation of the ΔS m and therefore also Δg requires precise determination of the melting conditions. The latter can be done using method of the liquid-solid coexistence, which use the fact that at melting conditions the solid and liquid phases remain in the thermodynamic equilibrium. Hence the special biphasic simulation box consisting 4096 RLM divided equally between separated by a small gap crystal and liquid phases. The box was created by the connection of the crystal structure registered during heating and the liquid box of the same high and wide. The long of the liquid part of the biphasic box was set to ensure density of liquid corresponding to that registered during cooling process. Subsequently, we simulated the biphasic system at various temperatures (keeping constant pressure) and determined the temperature at which the coexistence of both phases is observed for 10 ns (1 billion time-steps). The distinction between crystal and liquid was possible even by visual examination. When the system exhibits coexistence of two phases for 10 ns we recognize given temperature as the melting temperature (if needed we extended simulation time to complete melting or crystallization of the biphasic system). The determined melting temperatures equal 50 K, 36 K, 128 K, 65 K, 194 K for system A, B1, B2, C1, C2 respectively. According to Eq. (1), calculation of the nucleation rate requires also determination of the solid-liquid interfacial energy, which was estimated from the formula suggested in ref. 62 . At isobaric conditions of the reference pressure formula takes succeeding form , where typical 0 4 0 γ ≈ . was employed. The last quantity necessary to calculation of J is the diffusion constant, which was determined using GROMACS software from the means square displacement for long times. The obtained temperature dependence of D was subsequently approximated by the Vogel-Fulcher-Tammann equation.
The curves of the nucleation rates are presented in Fig. 3 for all studied systems. It is seen that the nucleation rate is the highest for the system C2, i.e., for the system which indeed crystalizes. Additionally, one can observed that around the temperature at which crystallization occurs ( = T K 130 , which is 64 K bellow the melting temperature) J for system C2 is 10 times higher than for C1 and B2 systems. This finding suggests that for other studied systems the nucleation rate is not sufficient for formation of the stable nuclei and then for the initiation of the crystallization process. Consequently, the crystallization does not take place for those systems. Moreover, we would like to note that the gain in the intermolecular attraction caused by increase in the dipole moment not always moves nucleation rate curves further from the melting temperature -maximum of J for system B1, which poses dipole moment (oriented alongside longer molecular axis), is closer to T m than the maximum of J for system A, i.e., for the system without dipole moment. This fact puts new insight on the straightforward connection between value of the dipole moment, related to it intermolecular attraction, and the crystallization tendency, which have been achieved from the studies of simple-liquids [19][20][21] . For these systems increase in the intermolecular attraction results in the shift of the nucleation rate further from T m . Consequently, the separation between curves of nucleation and crystal growth rates become higher, which impedes occurrence of the crystallization process. However, in the case of studied herein systems, one can observe in Fig. 3 that the nucleation rate curves are closer to the melting temperature and then to the respective crystal grow rate curves (crystal growth rate, U, were calculated in the same way like in ref. 20,21 , results not presented) for systems which do not crystalize. Consequently, the mismatch between nucleation and crystal grow, which is the highest for system C2 cannot be responsible for the occurrence of the crystallization process for the studied systems. On this base, one can suspect that the stability behavior of more complex systems than the simple liquids depends more significantly on the nucleation rate value than on the separation between J and U curves.
Moreover, comparing results obtained for all studied systems we can observe that playing the orientation and/or the value of the dipole moment it is possible to steer the nucleation rate. Creation of the dipole moment as well as the change of the existing one implies modification of the nucleation rate, see results for A, B1, and C1 presented in Fig. 3. However, it must be stressed that the obtainment of requested J behavior is a tricky task because molecular structure have to be considered as well. In should be also taken into account that too high the value of the dipole moment could lead to the overcoming of the previous decrease in J (caused by the given orientation of μ) and finally to the gain in nucleation rate (see results for B1 and B2 in Fig. 3).

conclusions
The simulations of the new model molecules, enable us to demonstrate that the interior orientation of the dipole moment has a significant role in the stability behavior of the liquid below the melting conditions. The change of the dipole moment arrangement could result in a total change of the tendency of the system to crystallization. Model molecules with the dipole moment oriented perpendicularly to their longest molecular axis exhibit higher tendency to the formulation of the crystal structure than those with the orientation of μ, which is parallel to the longest molecular axis. This result suggests that the impediment of the crystallization tendency, which is the effect of more complex molecular structure could be overcome by the interactions between properly arrange dipole moments. The determination of critical value of the dipole moment, which result in crystallization of RLM systems, requires further detailed studies. One can expect that this value depends on the interactions between atoms belonging to different molecules as well as on the ratio between diagonals of RLM molecules.
Interestingly, the presented herein results could be explained in term of the classical nucleation theory. The nucleation rate for the system with μ oriented perpendicularly the longest molecular axis is much higher than for system possessing μ of the same value but oriented parallel to the longest molecular axis. Consequently, the specific replacement of single atom within the molecule leading to the change of the dipole moment orientation could result in entirely different stability behavior of the system. Taking this fact into account, we can conclude that the molecular architecture, which is not considered in simulations of common model systems, has a crucial meaning for macroscopic properties of the system. At this point it must be noted that the structural anisotropy and the spatial charges distribution are not all components of the molecular architecture. Thus, the presented in this paper studies indicate on the completely unexplored direction for feature computational studies, which focus on the effect of given ingredient of molecular architecture in the structural, dynamic, and thermodynamic behaviors of real systems.