GPU-Accelerated Molecular Dynamics Simulation to Study Liquid Crystal Phase Transition Using Coarse-Grained Gay-Berne Anisotropic Potential

Gay-Berne (GB) potential is regarded as an accurate model in the simulation of anisotropic particles, especially for liquid crystal (LC) mesogens. However, its computational complexity leads to an extremely time-consuming process for large systems. Here, we developed a GPU-accelerated molecular dynamics (MD) simulation with coarse-grained GB potential implemented in GALAMOST package to investigate the LC phase transitions for mesogens in small molecules, main-chain or side-chain polymers. For identical mesogens in three different molecules, on cooling from fully isotropic melts, the small molecules form a single-domain smectic-B phase, while the main-chain LC polymers prefer a single-domain nematic phase as a result of connective restraints in neighboring mesogens. The phase transition of side-chain LC polymers undergoes a two-step process: nucleation of nematic islands and formation of multi-domain nematic texture. The particular behavior originates in the fact that the rotational orientation of the mesogenes is hindered by the polymer backbones. Both the global distribution and the local orientation of mesogens are critical for the phase transition of anisotropic particles. Furthermore, compared with the MD simulation in LAMMPS, our GPU-accelerated code is about 4 times faster than the GPU version of LAMMPS and at least 200 times faster than the CPU version of LAMMPS. This study clearly shows that GPU-accelerated MD simulation with GB potential in GALAMOST can efficiently handle systems with anisotropic particles and interactions, and accurately explore phase differences originated from molecular structures.


Introduction
Liquid crystals (LC) can be assembled from either small molecules [1,2] or polymers [3][4][5][6]. Compared with small molecular LC, liquid crystal polymers (LCP) offer both mesogen orientation and polymer mechanical merit [7,8]. The orientation of mesogens competing the penalty from the conformational entropy leads to abundant phase behaviors [9][10][11]. Nowadays, computer simulations have been widely used to improve our understanding of these phase behaviors which are not fully accessible by experimental approaches [12][13][14]. Adequate descriptions of the LC systems require an accurate model for both inter-and intra-molecular interactions, and also require system sizes that are large enough to observe the effects of molecular ordering [15]. A combination of coarse-grained Gay-Berne (GB) model [16][17][18] and graphics processing units (GPU)-accelerated algorithms [19][20][21][22] should be a highly promising way to provide simulation accuracy and efficiency in the study of larger systems considering more details of anisotropic particles, especially for LC phase transitions.
Coarse graining as one of the conceptual and technical ways that smoothes over or averages out some of fine details, has been widely used to simulate the longer time-and larger lengthscale dynamics and phase behaviors [23,24]. The coarse graining of LC through representing groups of atoms by single interaction sites leads to a remarkable acceleration in molecular simulation. Conventionally, LC mesogens are coarse grained as spheres in a wireframe graph, such as Kihara model [25], multisite model [13], and Lennard-Jones (LJ) model [26]. The GB potential employs soft ellipsoids rather than spherical particles to model particles with anisotropic characteristics [16,17], typically like discotic metallomesogens and dendrimers etc [27][28][29]. The interaction and the corresponding geometric packing of ellipsoidal particles have four forms, i.e., side-to-side, end-to-end, cross and T-shape. Through the consideration of anisotropic attraction and repulsion, the GB potential has been successfully used to model thermotropic LC phase transition for mesogens in small molecules [30][31][32][33][34]. De Miguel et al. identified the isotropic, nematic and smectic-B phases in small molecular LC phases and found that the stability of nematic phase is strongly influenced by the anisotropy in the well-depth of GB potential [35,36]. Luckhurst and his coworkers found that the stability of the isotropic and nematic phases is dominated by short range anisotropic repulsive forces, while the stability of the smectic-A phase is dominated by the anisotropic attractive forces in GB potential [37].
According to the location of mesogens in backbone or in side-chain pendant, LCP are classified into main-chain liquid crystal polymers (MCLCP) and side-chain liquid crystal polymers (SCLCP) [38]. The competition between the enthalpy gain from the global packing of mesogens and the entropy penalty from the conformation of polymer backbones and spacers leads to complicate phase behaviors which have attracted intensive studies [3,39,40]. Wilson et al. pioneered the GB/LJ model and offered a clear statement on the influence of odd-even spacer lengths on the LC phase transition in MCLCP [3,41,42]. In the GB/LJ model, rigid units are modeled by ellipsoidal particles described by GB potential, and flexible spacers are described by conventional LJ spherical particles. Allen et al. applied the GB/LJ model and proposed the mechanism for the thermotropic isotropic-LC phase transition of polyester with temperature decreasing [3]. Compared with well studied MCLCP, for SCLCP, the structure of polymer backbone, the grafting density, the pendant mesogens, and the flexible spacer length all exert important influence on the mesophase morphology [43]. SCLCP with highly flexible backbones (e.g., methacrylate based) and longer spacers have a lamellar smectic phase, while the same polymers with a short spacer typically have the nematic and isotropic phases at high temperatures [34]. Experimentally, neutron scattering [44], X-Ray Scattering [45] and birefringence measurements [46,47] have been steadily revealed abundant morphologies originated from the assembly of molecules with different mesogens and architectures. How to accurately address the relationship between the morphology of LC phase and the molecular structures is still of key importance.
Although GB potential has been wide accepted as an accurate model for LC mesogens, the computational complexity aroused challenges in simulations to identify LC phase transition in sufficient temporal and spatial scales. Various efforts including the simplification of energy function, the utilization of parallel computation or GPU-acceleration have been paid to reduce computational cost for GB potential computation. Zannoni and coworkers proposed a simplification by a sum of stretching and bending potentials for adjacent GB mesogens which can significantly reduce the computational cost [48,49]. A modified GB potential proposed by Persson can save 10% to 20% computational cost with reasonable approximation compared to the original GB potential [50]. Alternatively, Wilson et al. developed parallel MD simulation algorithms named GBMOL and GBMOL DD [39] to enlarge simulation systems. They realized a simulation with more 40960 atoms and observed novel patterns of LC phases. Meanwhile, the utilization of GPU acceleration enables MD simulation to study complicate phase behaviors of polymers within acceptable computational time [51][52][53][54]. Sunarso et al. used a GPUaccelerated MD simulation to study the generation of macroscopic flow through molecular reorientation in nematic LC phase which is triggered by an electric field [53]. They found that the GPU-acceleration can provide about 50 times faster than that using a single CPU core on an identical simulation. A series of software packages utilizing GPU-acceleration in generalized MD simulation have been developed, such as LAMMPS-GPU [55], HOOMD-blue [56] and GALAMOST [57]. Among them, GALAMOST fully took advantages of the mesoscopic simulation techniques including the soft anisotropic particle model [58], the iterative Boltzmann inversion numerical potential coarse-graining method [59], the hybrid particle-field MD [60], and the chain-growth polymerization model [61] etc. It has shown great potentials to efficiently model various nanoparticles [62][63][64] and polymeric materials [65,66].
In this study, we developed GPU-accelerated MD simulations equipped with coarse-grained GB potential implemented in GALAMOST package to investigate the LC phase transitions. The model and simulation details are presented in Section 2. Computation cost is compared in the simulation of mesogens in small molecules in Section 3. The difference in phase morphologies originated from mesogens in small molecular LC, MCLCP and SCLCP are discussed in Section 4.

Gay-Berne Potential
Generally, the interaction between spherical beads can be calculated by a LJ potential, U LJ where r ij is the distance between beads i and j located at r i and r j . The parameters ε 0 and σ 0 are taken as the units of energy and length, respectively. Then the GB potential is defined using an analogous form of the LJ potential The vector and parameter in the definition are depicted in Fig 1. Herer ij is the unit vector of r ij and û i is the principal vector of an ellipsoidal particle i. Functions εðû i ;û j ;r ij Þ and sðû i ;û j ;r ij Þ describe the interaction well-depth and the shape of two ellipsoids. If set one GB ellipsoid to present a typical calamitic mesogen, σ 0 is about 0.5Å and ε 0 is 1.38×10 -21 J. The shape function is defined as where the shape anisotropy parameter χ = (κ 2 -1)/(κ 2 +1) with the geometric aspect ratio κ = σ e / σ s . Here the σ e is the diameter of the major axis and the σ s (= σ 0 ) is the diameter of the minor axis of an ellipsoid. The interaction function εðû i ;û j ;r ij Þ is decoupled into and ε 0m ðû i ;û j ; where the interaction enhanced anisotropy parameter χ' is defined as w 0 ¼ ðk 01=m À 1Þ=ðk 01=m þ 1Þ. Here κ' = ε s / ε e represents the anisotropy of interactions and ε s and ε e are the well depths for the side-by-side and the end-to-end configurations, respectively. The empirical tuning exponents ν and μ are originally set to the values 1 and 2, respectively [3,16]. In order to reduce the computational complexity of GB potential, Persson introduced a modified GB (MGB) interaction [50]. Both the interaction and the shape functions are simply The interaction profiles of GB and MGB in four typical orientations are shown in Fig 1. The MGB has degenerate state profiles for side-by-side and cross interactions.

Coarse Grained Model
For the sake of a reliable comparison on the efficiency of GPU-accelerated MD simulation, we adopted a classical coarse grained model for mesogens proposed by de Miguel et al., where the temperature induced LC phase transition of small molecules has been clearly addressed [36].
The mesogens in small molecules, MCLCP and SCLCP are identical. It has a geometric aspect ratio κ = 3 and an interaction aspect ratio of κ' = 5. The potential is truncated at a distance r cut = (κ + 1)σ 0 and shifted the interaction at r cut to zero. An individual mesogen is treated as a linear rotor, and the moment of inertia around the principal axis of an ellipsoid is . This model has been used to yield the principal thermotropic LC phases [67]. For SCLCP, the main chain is connected by spherical beads (diameter of σ 0 ) and a flexible spacer. A harmonic interaction is applied to confine the length and the angle of a spacer while it does not occupy any volume. The parameters describing the GB/LJ interactions are obtained using the venerable Lorentz-Berthelot mixing rules [68]. We use a value σ e = 2.0 for the end-to-end manner, σ s = 1.0 for the side-by-side manner, ε s = 1.29 for the maximum well depth of the GB/LJ interaction and ε e = 0.58 for the minimum well depth of the GB/LJ interaction, which is the geometric mean of the GB and the LJ well depths [Fig 2(D)].
Based on this coarse-grained model, the energy function to guide MD simulation is written as where N LJ and N GB are number of LJ particles and GB particles, respectively. U ij LJ , U ij GB and U ij GB/LJ represent the nonbonded interaction potential. The harmonic interactions of bond and angle are used to model the connectivity between neighboring beads and the rigidity of polymer chains. N bond is the number of bonds, k 0 the spring constant, and r 0 the equilibrium bond length of flexible spacer, r is the bond length equal to the length of the site-site vector S ij . Chain rigidities are described by a bond bending potential between two consecutive bonds. N angle is the number of angles, k angle is the bending constant, θ is the bond angle and θ 0 the equilibrium bond angle. The spring constant k 0 = 100, the equilibrium bond length r 0 = 0.1, the bending constant k angle = 90 and θ 0 = 0.0. The backbone of SCLCP is considered as a bead-spring chain with k 0 = 100 and r 0 = 1.0, where r 0 is the equilibrium bond length between the center of mass of two adjacent beads.

Simulation Settings
To study the LC phase transition of mesogens in different molecules, simulations are in NVT ensemble with periodic boundary condition is used. For mesogens in small molecular LC systems, the simulations are performed with 1000 mesogens. For MCLCP systems, they contain 100 polymer chains and each has 10 mesogens. For SCLCP systems, there are 10 mesogens and 19 beads in 64 polymer chains. These three system have the same volume (V = 15σ 0 ×15σ 0 ×15σ 0 ) and volume fraction of particles (ρ Ã %0.47). Further, to confirm finite size effect in the simulation systems, we also carried out simulation in the simulation box with double size. The large systems contain 8000 mesogens for small molecular LC and MCLCP with 20 mesogens in each chain, as well as 4960 mesogens and 9672 beads for SCLCP with 20 mesogens and 39 beads in each chain. Ellipsoids are mesogens, spheres are backbone connectors in SCLCP and the springs demonstrate harmonic interactions. θ is the angle between the major axis of the two adjacent GB mesogens (u i , u j ) for MCLCP (C), while it is the angle between the major axis of the GB mesogen (u i ) and the vector from the center of the GB mesogens to the adjacent LJ beads on backbone (r ij ) for SCLCP (D). The site-site vector S ij connects the two adjacent GB/GB sites or GB/ LJ sites placed in terminal position for MCLCP and SCLCP, respectively.
A standard leap-frog algorithm for anisotropic systems is used to solve the equations of motion, and the dimensionless MD time step dt Ã = (mσ 0 2 /ε 0 ) 1/2 dt = 0.001. The reduced temperature T Ã = k B T/ε 0 is controlled using Nose-Hoover thermostat with a relaxation time τ Τ = 0.1, where k B is the Boltzmann constant. To observe the phase behavior of the model systems, the simulated annealing process is applied to gradually cool down the systems from equilibrated isotropic melts. Systems are equilibrated for at least 1×10 7 and 2×10 7 time steps for mesogens in small molecules and polymers, respectively. The mean properties are represented by statistical averages over ten samples evenly taken from the last 3×10 6 time steps. To test the simulation cost in small molecule LC systems, the number of particles range from 1×10 3 to 1×10 6 , and accordingly system sizes from 10σ 0 ×10σ 0 ×30σ 0 to 100σ 0 ×100σ 0 ×300σ 0 to keep the constant number density ρ = N/V = 0.33. All simulations using either GALAMOST or LAMMPS are performed in singleprecision.
Simulations were carried out using the GPU-accelerated package GALAMOST [57] and the well known LAMMPS package [69] for simulation cost comparison. A server equipped four NVIDIA Tesla K20C GPUs (each 2496 cores) and two Intel Xeon E5-2687w CPUs with 128G RAM is used. Each simulation trajectory utilizes either one GPU or CPU.

Characteristic Parameters
The global orientation of mesogens in the simulation system are viewed from the order parameter S, which is computed by 3 The value of S is in the range of 0 to 1, 0 for random distributed and 1 for completely orientated mesogens. The angle γ i is between the major axis of the mesogen i and the principal axis of all mesogens. The principal axis is determined from the largest positive eigenvalue of the secondrank symmetric tensor Q 3,41 , where the unit vector u i points along the major axis of molecule i. The Kronecker delta function δ αβ is 1 when α and β are identical, and 0 otherwise. The spatial location of mesogens are indicated from the radial distribution functions g(r) of mesogens, which gives the probability of finding a particle in the distance r from another particle 41 , The range of g(r) is from 0 to large positive number. It equals 1 indicates random distribution, larger value for enrichment and less than 1 for depletion of mesogens at given separation range.
Then the orientational correlation of mesogens at given separation is presented by g 2 (r), defined as 42 The value of g 2 (r) equals to 1 corresponds to that all mesogens are in side-by-side or end-toend packing and a value close to 0 indicates mesogens are in isotropic phase. The association of mesogens in LC systems can be viewed by the second virial coefficient A 2 , which is calculated using the definition of [70] where N ij and U(i, j) is the total number of interaction pairs and the interaction energy between two particles i and j, respectively. A negative A 2 indicates that mesogens prefer close packing, and a positive one corresponds to dispersed state.

Protocol and Implementation in GALAMOST
To eliminate the time-consuming data transfer between host memory and device memory, GALAMOST is designed to implement all MD computations on GPUs, including integration, building neighbor list, and the calculation of non-bonded and bonded interactions, etc. By continuous optimization, GALAMOST can run efficiently on a single GPU. Some advanced optimization techniques are employed, such as the SFCPACK sorting method which rearranges the order of data of particles stored in host and device memories to increase the probability of cache hit at memory reading. In addition, the GPU-accelerated neighbor list building algorithm is used to further improve the computational efficiency. We employ Compute Unified Device Architecture (CUDA) to harness the computational power of GPUs. All GPU threads execute the same instruction, while process different data. The execution mode is thereby called single-instruction-multiple-data (SIMD). We usually design the GPU-accelerated MD algorithm by the principle that one GPU thread tackles the computational task of one particle. GB interactions as non-bonded interactions are calculated in "one thread one particle" mode. The algorithm is to use a neighbor list that lists all interacting GB particles for each particle, built beforehand. Because of the independence of parallel GPU threads, a pair of interacting particles is inevitably separated in neighbor list. The GB forces and torques for each particle are calculated and summed by a corresponding thread on the GPU via searching the particles within cutoff in neighbor list process.

Compare Simulation Cost and Accuracy of MD Simulation on Small Molecular Liquid Crystals
Simulation cost is compared between LAMMPS and GALAMOST in the original and modified GB potential forms, with systems spanning from 10 3 to 10 6 particles. The average simulation time over 10 parallel runs against the number of mesogens is presented in Fig 3. For small system (N<10 3 ), the GPU-acceleration is not prominent compared with conventional CPU-based algorithm. When the number of particles reaches 10 6 , the cost of GPU MD simulation using GALAMOST were about 4 times and more than 200 times less compared the GPU and CPU MD simulation in LAMMPS, respectively. Compared with the original GB potential, the MGB potential which takes a simplified functional form affords about 30% to 50% in the reduction of simulation cost.
We also decoupled the three major time-consuming processes and their comparisons for a simulation system with 10 6 particles are shown in Fig 4. The GPU-accelerated simulation show advantages in all the three functions of force, neighbor list and integration processes. Because of the high complexity of GB potential, force calculation turns to be the major time consuming process. The success of GPU acceleration may benefit from two important aspects: (1) the computational power of a single GPU, such as Tesla K20C with a theoretical peak 3.95 Teraflops of single precision computation throughput, is usually more than a hundred times faster than that of a CPU core, such as Xeon E5-2687w with 21.56 Gigaflops per core; (2) GALA-MOST is highly optimized such as the memory sorting technique which can significantly reduce the time for nonbonded force calculation, and the elimination of inter-memory communication [57].
The snapshots of simulation configurations under equilibrium from low to high temperatures, associated with the phase transition of small molecular LC are presented in Fig 5. At low temperature (T Ã = 0.6), mesogens are condensed in a smectic-B phase. The global orientation of mesogens is conservative for the simulation guided by GB potential, independent on either the CPU-approach used in LAMMPS or the GPU-acceleration approach implemented in GALAMOST. At high temperature (T Ã = 0.8), there is no global orientation, exhibiting typical isotropic phase. The configurations based on GB interaction produced by the two simulation packages are non-distinguishable, while the simulation guided by MGB which has degenerate states of side-by-side and cross interactions did not show the orientational transition at this temperature window.
Further comparison of simulation accuracy through the global order parameter S is shown in Fig 6. The temperature dependence of the order parameter from the configurations generated by GB potential with GALAMOST and LAMMPS packages are almost identical. At low temperature, S levels off at 0.96, indicating that measogens are highly orientated. At a critical temperature at T Ã = 0.66, there is a smectic-isotropic phase transition. At high temperature (T Ã >0.9), the value of S is close to 0, meaning that globular orientation of mesogens vanishes. However, the phase transition guided by MGB spanned in a much broader temperature window. It suggests that the MGB potential is not suitable for LC phase transition study although it has merits in simulation cost saving.
To further confirm the accuracy of GALAMOST and GB potential in the simulation of LC phase transitions, we reproduced the phase diagram (Fig 7) using the identical simulation systems reported by de Miguel et al. [36]. The data points to identify the phase boundary are determined from configurations generated in simulated annealing processes. Either the global order parameter S sharply increases or the g(r) profile has new peaks (it is also combined with direct visualization of configurations) is regarded as a phase transition point. Using this criterion, the GPU-accelerated MD simulation equipped with GB potential can reliably and accurately repeat the phase diagram. Therefore, in the following section, the LC phase transitions for mesogens in different molecules are only simulated by the GALAMOST package with GB potential.

Liquid Crystal Phase Transition of Mesogens in Different Molecules
Typical snapshots of simulation configurations for small molecules, MCLCP and SCLCP at ordered, intermediate and disordered states with temperature increasing are shown in Fig 8. At high temperature, the location and the orientational orders disappeared, all systems exhibit typically isotropic phases. At low temperature, mesogens in small molecules are packed in smectic phases with layered and ordered structures, while nematic phases appear in both SCLCP and MCLCP systems. Although small molecules and MCLCP have the same number of mesogens, the competition between the orientational ordering of mesogens and the entropy maximization of the backbone in MCLCP leads to higher transition temperature than that in small molecular systems. Compared with the one-step transition from liquid to ordered LC phase for mesogens in small molecules and MCLCP, SCLCP undergo two-step for this transition as temperature decreases. When temperature decreases from 1.0 to 0.7, mesogens in SCLCP gradually assembled into enriched domains. In this process, nematic nuclei are formed and rapidly grow in size (see S1 Fig). Further decreasing to low temperature (T Ã = 0.4), SCLCP forms multi-domain nematic phase where mesogens in each domain has a favored orientation. It is due to that the backbones and spacers hinder the rotational orientation of the mesogenes. The phase transition for SCLCP in this simulation is consistent with the results in birefringence and optical transmission experiments [46,47]. Boamfa et al. studied various SCLCPs and found that during the first-order isotropic-nematic phase transition, the favored orientation for each domain is determined by the ordering fluctuation which originates from the polymer backbone coupling effect.
The orientational order parameter S for mesogens in small molecules, SCLCP and MCLCP systems are shown in Fig 9(a) and 9(b). The temperature dependence of S for both small molecules and MCLCP follow a Sigmoidal curve [71] associated with the LC phase transition. Compared to the sharp transition of small molecular LC, MCLCP spans a broader temperature window and the critical transition temperature is at a higher temperature T Ã = 12.0. The LC transition of MCLCP may experience an coexistence of anisotropic and isotropic phases, as experimentally revealed by Shilov et al. [72] The higher critical temperature for MCLCP results from the restricted translational and rotational motions of the mesogens by the rigid backbone, and the additional energy penalty derived from the orientation of polymer backbones [73]. SCLCP shows multi-domain morphology with local orientation as indicated from the global order parameter S is nearly 0. Correspondingly, in experiment, the birefringence of this multidomain morphology is very small due to the absence of a common director orientation. Furthermore, in order to show the degree of local orientation, we propose the probability P to count the preference of local ordering composed by neighboring mesogens. The probability P is defined as proportion of the angle between their major axes satisfies |cos(u i Áu j )|>0.8 outcomes to the total number when two mesogens with separation distance less than 3.4σ 0 . For small molecules and MCLCP systems, the temperature dependence of P matches S well. However, as SCLCP only has local orientation, S does not reveal the corresponding ordering associated with the LC transition, while P clearly shows such transition.
We also calculated the second virial coefficient A 2 to view the association preference of mesogens, as shown in Fig 9(c) and 9(d). Increasing temperature resulting in the decrease of -A 2 , which can be seen in all three systems. It may be attributed to the decrease of attractive interaction and the reduction of closely packed neighboring molecules according to the definition of A 2 . Further, due to the extended conformation of rigid main-chain and higher temperature zone, A 2 of MCLCP is much smaller than those in small molecular and SCLCP systems. At low temperatures, strong interactions among mesogens can assemble almost all molecules together, and mesogens in all three systems show an ordered structure and a close-packed arrangement. At high temperatures, the asymptotic value of A 2 is close to 0 in small molecular and SCLCP systems, suggesting the ordered packing of mesogens is not favored. In other words, mesogens are in isotropic state. However, the A 2 in MCLCP systems becomes positive. This is due to that the location of mesogens in MCLCP system is disordered, and mesogens prefer to orient along the rigid main chain and repel each other. It has also been reported by Withers et al. [74] that the second virial coefficient gradually increases to positive with temperature increasing for rigid polymers. Increasing temperature weakens the attractive interaction associated with the packing of orientated mesogens, and the rigidity of the MCLCP backbone leads to slightly repulsion among mesogens.
To further clarify the distribution and the orientation of mesogens, the radial distribution function g(r) and the orientational correlation function g 2 (r) are presented in Fig 10. All systems show an obvious first correlation peak at 1.12σ 0 which corresponds to side-by-side and cross packing of mesogens, and the former is the majority as shown from the positive peak values of g 2 (r). The second and third peaks are two and three times the separation distance of the first peak, respectively. It indicates that the packing of mesogens achieve long-range ordering. Further decreasing temperature makes the correlation peaks more prominent, in good agreement with that lower temperature leads to higher ordering of thermotropic LC. The value of g 2 (r) at long distances tends to S 2 , which is consistent with the results of simulation and theory 3 . Different from mesogens in small molecules, there is no sign of long-range positional order for MCLCP at low temperature (T Ã = 7.0). The curve of g 2 (r) shows the nearest-neighbor peak (r/σ 0 %1.12) and level-off at S 2 when mesogens are separated far enough. It indicated that the system appears to form a nematic, not a smectic phase. In the isotropic phase (T Ã = 13.0), g 2 (r) decays to zero at large separations and a small peak (r/σ 0 %3) reflects the strong intramolecular correlations. Compared with uniform orientation in nematic phase, there is no sign of longrange orientational order. For SCLCP, at low temperature (T Ã = 0.4), the curve of g(r) show a peak at r/σ 0 %1.12 for the closely packed neighboring mesogens. There is no obvious peak at 3σ 0 because the end-to-end packing of mesogens is inhibited by the backbone in SCLCP. Compared both g(r) and g 2 (r) for mesogens in polymers with those in small molecules, the connection of polymer chain can efficiently screen the long range correlation of orientated mesogens. Finally, in LC phase, the side-by-side and the end-to-end are favorable, and the cross and the T-shape have low occurrences.
Additionally, the finite size effect on the simulation of LC phase transition is also studied. The results from the simulation of large systems are presented in S2 and S3 Figs. Compared with the results shown above, it can be seen that enlarging simulation box only slightly changes the absolute values. The order parameter S, the probability of the local orientation P, the secondary Virial coefficients A 2 and the correlation functions all show similar trend. These results suggest that the simulation system in this work is large enough to present reliable description for LC phase transition for mesogens in different molecules. Meanwhile, we also found that longer polymer chain causes MCLCP to form ordered structure at lower temperature, while SCLCP has phase transition at slightly higher temperature.
In order to give a deep understanding in phase transition in LCP systems, the change of conformations for individual polymer chains along with temperature, characterized by the average end-to-end distance (R f ) for polymer backbones in MCLCP and SCLCP systems, is presented in Fig 11. The R f of MCLCP sharply decreases from T Ã = 12.0 to T Ã = 13.0, suggesting that chain collapse is coupled with the isotropic-nematic phase transition. The chain collapse results in less orientation among intra-chain mesogens, as revealed from the decrease in g 2 (r). In contrast, on cooling, the shape of the backbone of SCLCP is always a coil with R f only slightly increasing, which is associated with the multi-domain isotropic-nematic transition. Such increasing originates from that polymer chains spanned multiple orientated domains. In additional, assuming that polymers are Gaussian, the ideal R f values are 9.5σ 0 and 4.4σ 0 for MCLCP and SCLCP, respectively. The actual R f values under all conditions are significantly larger, indicating that the presence of mesogens extends polymer backbones.

Conclusion
We developed GPU-accelerated Molecular Dynamics simulations combined with coarsegrained Gay-Berne potential to investigate the temperature induced liquid crystal phase transition of mesogens in small molecules, main-chain liquid crystal polymers and side-chain liquid crystal polymers. On cooling, mesogens in small molecular liquid crystals exhibit the isotropicsmectic transition, main-chain liquid crystal polymers assemble into nematic phases, and mesogens in side-chain liquid crystal polymers undergo a two-step process including nucleation of nematic islands and formation of multi-domain nematic texture. Mesogens in small molecules and side-chain liquid crystal polymers share lower phase transition temperature, and main-chain liquid crystal polymers need much higher temperature to overcome the energy barrier from the orientation of rigid polymer backbones.
For systems containing up to one million particles, the GPU-accelerated simulation in GALAMOST is about 4 times less computational cost than the GPU simulation, and at least 200 times faster than the conventional CPU simulation in LAMMPS. The acceleration is The radial distribution function g(r) and the orientational correlation functions g 2 (r) at various temperatures for small molecular, MCLCP and SCLCP systems. The three vertical dash lines label the location of the side-by-side and cross, the 2 nd nearest side-by-side and T-shape, and the end-to-end packing of mesogens from left to right.  mainly contributed from the efficient handling of the high computation complexity of Gay-Berne potential, which is the major time-consuming process. The modified form of Gay-Berne potential is confirmed not suitable in study of liquid crystal phase transition though it has merit to reduce the computational complexity. The GPU-accelerated MD simulation with GB potential implemented in GALAMOST package is efficient and accurate to study the complex phase transition in liquid crystal systems. It provides a useful tool to simulate anisotropic particles and interactions for larger system size over longer time.