Development of Coarse-Grained Liquid-Crystal Polymer Model with Efficient Electrostatic Interaction: Toward Molecular Dynamics Simulations of Electroactive Materials.

Liquid-crystal polymers (LCPs) are well known materials for functional sensor and actuators, because of their high-responsiveness to an electric field. Owing to their complex physical nature, however, the prediction of the functions of LCPs is a challenge. To attack this problem from a molecular point of view, a simulation study is a promising approach. In this work, for future applications of molecular dynamics simulations to problems involving an electric field, we develop an LCP model which consists of coarse-grained mesogenic molecules and smeared charges. For the smearing function of the electrostatic force, the Gauss error function is introduced. This smearing is optimized to attain a reasonable accuracy for phase transition phenomena of liquid crystal while numerical instabilities arising from the singularity of the Coulomb potential are circumvented. For swelling systems, our LCP model exhibits the characteristics of both liquid crystals and unentangled polymer chains; orientational order of the mesogenic units and Rouse-like relaxation dynamics. Our coarse-grained LCP model successfully incorporates electric charges and dipoles and is therefore applicable to problems concerning an electric field.

1. Introduction

Liquid-crystal polymers (LCPs) consist of open chain compounds containing mesogenic units. It is well known that LCPs exhibiting liquid crystal phases show high-responsiveness to external stimuli such as an electric and magnetic fields, and irradiated light [1,2,3,4]. Cross-linked LCPs, referred to as liquid crystal elastomers (LCEs), liquid crystal networks (LCNs) and so on, have elastic characteristics, while maintaining the above-mentioned responsiveness that, in combination with cross-linked network, gives rise to macroscopic deformations [5,6,7,8,9,10]. This unique property of cross-linked LCPs indicates their applicability to sensing devices as well as actuator devices [11,12,13]. The development of these devices, however, is a challenge because the molecular-level mechanics of the deformation is not understood enough. The complex structures, dynamics, and physical properties of LCPs at different time and length scales over a wide range are closely related. In other words, a slight modification of atomistic-scale properties such as molecular architectures can drastically change meso- and macro-scale characteristics of LCPs; the understanding of the behavior of LCPs is a typical multiscale problem. Therefore, to attack this critical problem and predict physical properties of LCPs, simulational studies from molecular point of view should be introduced. Molecular models and simulations for LCPs have been developed and performed over ten years [14,15,16,17].

Skačej et al., developed an off-lattice coarse-grained LCP model using the soft-core Gay-Berne (SCGB) ellipsoidal potential, and carried out Monte-Carlo simulations under simplified body-force-like electric field conditions [16]. They observed the reorientation of ellipsoids induced by the external field, and the resulting deformation of LCPs. Whitmer et al., developed a flexibly bonded Gay-Berne (GB) model, and carried out molecular dynamics (MD) simulations [17]. A polydomain-monodomain transition of LCPs was observed under uniaxial strain. These previous studies indicate that coarse-grained simulations can provide a promising tool for the investigation of mesoscopic deformations of LCPs induced by external stimuli. It is important to note, however, that molecular simulations for macroscopic deformation of LCPs are still beyond the capability of current computational resources. Thus, hierarchical linkage between molecular simulations and macroscopic methods, such as the finite element method (FEM), is required. In FEM, electric-field-induced deformations can be computed using the piezoelectric tensor analysis [18]. For the linkage between two simulations with different time and length scales, therefore, the piezoelectric tensor should be directly calculated from molecular simulations. To enable the direct calculation, the exact expression of the dielectric character of molecular systems is necessary. Placing point charges on the coordinate of atoms is one simple and efficient way to express the dielectric character of molecular systems, and is a normal way for all-atomistic molecular simulations. For coarse-grained models, however, placing point charges has difficulty because of the bad combination between softened interaction for coarse-grained models and Coulomb interaction for point charges. Namely, the Coulomb potential has a singularity at zero distance of interaction while the softened potential does not. To circumvent the problem, a smearing function which removes the singularity can be introduced for the Coulomb potential. This “smeared charge” has been used for several simple coarse-grained models [19,20], but has not been developed for complex coarse-grained LCP models.

In this work, toward future molecular simulations of LCPs, we develop a coarse-grained LCP model. To enable the direct calculation of piezoelectric tensor, the dipole moment on a mesogenic unit is expressed using a pair of point charges. The Coulomb force for point charges is smeared using the Gauss error function. With using this smeared charge, the coarse-grained molecular simulations are stabilized while maintaining reasonable accuracy to investigate the phase transition phenomena of LCPs qualitatively. We perform MD simulations using this model for LCP swelling systems without the effect of an external field. The results show that our model exhibits the characteristics of liquid-crystalline orientational order and Rouse-like relaxation dynamics of unentangled polymer chains. We conclude that our coarse-grained LCP model has the possibility of hierarchically linking molecular and macroscale simulations.

2. Methodology

2.1. Coarse-Grained LCP Model

Our coarse-grained LCP model can be thought of as a modified version of the model of Skačej et al. [16,21]. The main difference from their original model is the introduction of point charges for the direct calculation of the piezoelectric tensor. A mesogenic dipole moment is represented by a pair of charges placed on the SCGB ellipsoid. For stabilizing coarse-grained molecular simulations, the electrostatic interaction between point charges should be carefully treated. The Coulomb potential cannot be used as is together with soft-core potentials, because the former has a singularity at zero distance of interaction while the latter does not. Without a proper treatment for point charges, simulations become unstable. The introduction of a smeared charge is one of the possible ways to handle this problem. Several types of smeared charges have been introduced for dissipative particle dynamics simulations [19,20]. In this work, the Coulomb force between point charges Q and q is modified by introducing the Gauss error function in the following manner: where is the electrostatic potential, k is Coulomb constant that will be set to unity in our simulations, is the vector representing the relative distance between two charges, , and in the error function concerns the strength of smearing. The following four conditions are desirable for the smearing function. (i) The function is continuous at any r. (ii) The singularity of Coulomb force at can be removed. (iii) The shape of the original Coulomb force is kept, except near . (iv) The implementation and computation are easy. Since the Gauss error function with large satisfies the above four conditions, it is ideal for a modified Coulomb force and useful for a wide variety of coarse-grained models.

The implementation for mesogenic ellipsoid follows that of the course-grained MD program COGNAC [22,23]. As shown in Figure 1, a GB ellipsoid is represented by two segments that are placed at the ends of its principal axis. For simplifying the computation, a pair of point charges is placed on the two segments belonging to the same ellipsoidal unit. The segments are bonded so that they constitute straight chain polymers. COGNAC is developed as a general purpose coarse-grained molecular dynamics package to adopt various kinds of coarse-grained model and to be able to introduce the user-defined potential function. Thus, the coarse-grained model used in this study can be easily implemented.

Figure 1

Schematic illustration of our coarse-grained liquid-crystal polymer (LCP) model.

For the intermolecular interaction between ellipsoids, SCGB potential [16] is used. The original GB potential is a well-known anisotropic interaction potential which has been widely used to study various kinds of liquid crystal phases and systems [24]. The GB potential energy is expressed as follows: where represents the vector between the centers of the mass for a pair of ellipsoids, i and j, and . The unit vector is along the long axis of the ellipsoid i. The length parameter for the ellipsoid pair is given by where

Here, , with and denoting the range parameter for the end-to-end and side-by-side configurations, respectively. The anisotropic energy parameter is written as

Here, is a characteristic GB interaction strength, and are dimensionless parameters, and where , with and denoting the potential well depth in the side-by-side and end-to-end configurations, respectively, and is set to as a characteristic van der Waals diameter.

In the SCGB potential, the strong core repulsion in the original GB potential is replaced by a weaker linear repulsion. The SC potential energy is written as

Here, m stands for the potential slope and and are already defined above. The GB and SC potentials are merged using a sigmoidal logistic function so that the resulting potential reads

In this study, the van der Waals diameter , GB interaction strength and mass of ellipsoids are set to unity, and is set equal to . The potential parameters used are , , , , , . The cutoff distance is .

For the bonding potential between adjacent segments k and , a harmonic potential is introduced as below:

Here, k is the spring constant, r is the bond length and is the equilibrium bond length. This simulation system possesses two kinds of bonding, inter- and intra-ellipsoid bonding, since the GB ellipsoid is defined by two segments. For inter-ellipsoid bonding, and are adopted. On the other hand, for intra-ellipsoid bonding, and . We choose the spring constant for the intra-ellipsoid bonding much larger than that for the inter-ellipsoid bonding so as to prevent the deviation of the segment location from the end of the ellipsoids.

2.2. Simulation Detail

Our MD simulation system is a cubic cell of volume V including 121 LCP molecules each made up of 30 GB particles, and 3630 GB monomer particles. Therefore GB particle are in our system in total, and the number of GB particles constituting LCP molecules is equal to that of GB monomer particles. The dimensionless density is . The equations of motion are solved by the velocity Verlet algorithm [25] with a time step of . In the simulations, the NVE ensemble is employed and velocities of mass points are scaled every 100 time steps to adjust temperature. The simulation is started with randomly oriented configuration at high temperature , and the temperature is decreased stepwise with the interval . At each temperature step, at least time steps of the calculation are carried out so as to obtain relaxed states. The relaxed configuration is utilized as the initial structure of the next temperature step. For the evaluation of the relaxation time, at least time steps of calculations are carried out, and in addition, 5 samples with different initial configurations are averaged. All of the simulations are done by COGNAC [22].

3. Results and Discussions

The smearing strength of the electrostatic potential is adjustable by the parameter in Equation (1). Figure 2 shows the potential profiles for various values of , where . The original unsmeared potential is denoted by a dashed line. The potential is smeared more strongly with decreasing , but the effect of the smearing should be as small as possible due to its artificiality. On the other hand, a certain level of smearing is essential for the stability of the simulations. In order to find an appropriate value of , we consider the influence of on the phase transition behaviors. The phase transition behavior is evaluated by the temperature dependence of the order parameter S defined as where represents the number of GB particles composing the LCPs (not including the solvent GB particles), and the director denotes the average orientation direction of GB particles composing the LCPs.

Figure 2

The potential profiles for the smeared electrostatic potentials with the strengths . The original potential is given by a dashed line.

Figure 3 shows the temperature dependence of the order parameter at various smearing strengths, , for . For , the calculation collapses. For the system without electric charges (dashed line), the drastic temperature variation of the order parameter (nematic-isotropic phase transition) is observed around . With increasing , the phase transition behaviors are smoothed and the transition temperatures are shifted higher. The two plots for and are, however, completely overlapped, which means that the introduced error function does not affect the phase transition behavior at . In this region of , the potential profiles are indistinguishable from the original potential at as shown in Figure 2. From these facts, we conclude that the smearing of the electrostatic potential at does not affect the phase transition behavior, whereas it stabilizes the calculation. In other words, the potential profile should be set to at for preventing artificial effects on the relaxation state, and smeared enough at for stable calculations. Therefore, we adopt in the simulations below.

Figure 3

Temperature dependences of the order parameter for the smeared electrostatic potentials with the strengths . The order parameter without the electrostatic potential is given by a dashed line.

There are two possible reasons for the smoothed phase transition behavior by introducing the electrostatic potentials. One is the influence of finite system size. The phase transition behaviors are confirmed to be little affected when the linear dimension of the system is doubled, and therefore the finite size effect is not likely to be the major cause of the smoothed phase transition behavior in the presence of electrostatic interactions. The other is the alignment of ellipsoids induced by the electrostatic potential. The ellipsoids with large spontaneous polarizations are stabilized strongly due to the head-to-tail dipolar interactions that affect the phase transition temperature [26,27]. Moreover, when an electric field is applied throughout the nematic-isotropic phase transition, the temperature dependence of the order parameter is smoothed and the (smoothed) transition is shifted to higher temperature due to the electric-field-induced stabilization of the liquid crystal alignment [28]. In our model, the alignment of ellipsoids with strong polarization also generates global polarization, which is a situation similar to that of a system polarized by external electric fields. By introducing the electrostatic interaction, the nematic phase is stabilized especially at , where the isotropic phase appears in the absence of electrostatic interactions. As a result, the phase transition is shifted to higher temperature and smoothed by introducing the electrostatic charges.

Figure 4 shows the phase transition behaviors for several values of the electric charges, , where . For , the calculation is extremely unstable. With increasing the charge, the phase transition is smoothed and shifted to a higher temperature due to the electrostatic stabilization of the ellipsoid alignment stated above.

Figure 4

Temperature dependence of the order parameter for , where .

The understanding of the relaxation dynamics is crucial for the elucidation of the mechanical characteristics of LCPs. For this purpose, we analyze the autocorrelation of normal modes at the isotropic phase (), which is written as [29], with where p is the index of the normal mode, is the relaxation time of and represents the position of the center of the mass in the i-th ellipsoid. Figure 5 shows the mode index p dependence of the relaxation time for and . We find that the relaxation time do not depend on the strength of the electrostatic interaction in the isotropic phase. Interestingly, the relaxation of LCPs exhibits the Rouse behavior, where is proportional to [30]. This result suggests that the LCPs in our coarse-grained model are unentangled and flexible enough in spite of the rigid nature of the GB ellipsoids. This is partly because the swelling monomers improve the mobility of LCP chains and disturb the inter-molecular interaction. Furthermore, our LCP chains are likely to be too short to induce entanglement; more than around 100 monomers per chain are needed for the entanglement in the case of the typical bead-spring model [29]. Therefore, we speculate that the entanglement is observable in our LCP model by decreasing the solvent concentration and preparing longer chains.

Figure 5

The mode index p dependence of the relaxation time for the autocorrelation function of the normal mode upon the isotropic state, .

4. Conclusions

We developed an LCP model which uses coarse-grained molecules and smeared electrostatic charge toward the prediction of electrostrictive dynamics. For the smearing function of the electrostatic force, we introduced the Gauss error function. The smearing strength was optimized so as to attain a reasonable accuracy for phase transition phenomena of liquid crystals. The nematic-isotropic phase transition was smoothed and shifted to higher temperature due to the electrostatic interaction. This is because the strong polarizations on the aligned ellipsoids generate a global polarization, which stabilize the nematic phase. In addition, the relaxation analysis for the autocorrelation of normal modes revealed that our model exhibited unentangled polymer characteristics. In other words, we confirmed that our LCP model with smeared charges exhibited both characteristics of liquid crystals and polymers. Thus, we succeeded in the development of an LCP model that can react to the electric field, and our model will be useful for predicting the mechanical properties of LCPs induced by the alignment of the mesogenic units. In further studies, the piezoelectric tensor will be calculated by this model, which will enable the hierarchical linkage between microscopic modeling and macroscopic electrostrictive simulations. We believe that this linkage for the prediction of the macroscopic functions of LCPs leads to successful molecular designing for the development of LCP sensors and actuators.