Transferability of Coarse-Grained Force Field for nCB Liquid Crystal

Article pubs.acs.org/JPCB

Transferability of Coarse-Grained Force Field for nCB Liquid Crystal Systems Jianguo Zhang and Hongxia Guo* Beijing National Laboratory for Molecular Sciences, State Key Laboratory of Polymer Physics and Chemistry, Institute of Chemistry, Chinese Academy of Sciences, Beijing 100190, China S Supporting Information *

ABSTRACT: In this paper, the transferability of the coarsegrained (CG) force field originally developed for the liquid crystal (LC) molecule 5CB (Zhang et al. J. Phys. Chem. B 2012, 116, 2075−2089) was investigated by its homologues 6CB and 8CB molecules. Note that, to construct the 5CB CG force field, we combined the structure-based and thermodynamic quantities-based methods and at the same time attempted to use several fragment molecular systems to derive the CG nonbonded interaction parameters. The resultant 5CB CG force field exhibits a good transferability to some extent. For example, not only the experimental densities, the local packing of atom groups, and the antiparallel arrangements of nearest neighboring molecules, but also the unique LC mesophases as well as the nematic−isotropic phase transition temperatures of 6CB and 8CB were reproduced. Meanwhile, the limitations of this 5CB CG force field were also observed, such as the phase transition from nematic to smectic was postponed to the lower temperature and the resulting smectic phase structure is singlelayer-like instead of partially interdigitated bilayer-like as observed in underlying atomistic model. Apparently, more attention should be paid when applying a CG force field to the state point which is quite different from which the force field is explicitly parametrized for. The origin of the above limitations can be potentially traced back to the inherent simplifications and some approximations often adopted in the creation process of CG force field, for example, choosing symmetric CG potentials which do not explicitly include electrostatic interactions and are parametrized by reproducing the target properties of the specific nematic 5CB phase at 300 K and 1 atm, as well as using soft nonbonded potential and excluding torsion barriers. Moreover, although by construction this CG force field could inevitably incorporate both thermodynamic and local structural information on the nematic 5CB phase, the anisotropic diffusion coefficient ratios for different LC phases in both 6CB and 8CB systems are reproduced well. All these findings suggest that the multiproperty parametrization route together with fragment-based method provides a new approach to maximize the possibility to simultaneously reproduce multiple physical properties of a given molecule or related molecules with similar chemical structures at other state points.

1. INTRODUCTION The coarse-grained (CG) models can easily access the processes taking place at the mesoscopic level, such as the phase behavior and self-assembly, as the reduction of degrees of freedom of molecules in the CG models results in a considerable speed up (by a factor of at least 101) when comparing with the corresponding atomistic models. In consequence, how to develop the CG models which retain close connections to the underlying atomistic representation and can be used to make quantitative predictions for specific materials has attracted significant attention in recent years. In general, one should construct different CG models according to the target properties wanted to reproduce. Moreover, to reproduce the target properties accurately, a reasonable mapping scheme and an optimal CG force field (FF) are needed. The mapping of the atomistic structure to the CG model is not unique but at the same time very crucial in terms of both predicting power and efficiency. However, in the previous © 2014 American Chemical Society

literature of CG simulations, the procedure of determining suitable force field parameters has been emphasized, and few studies have pursued with the aim to design rational mapping schemes and learn something about their capacities. As addressed in our previous work,1 to reproduce the phase behavior of the 5CB (4-cyano-4′-pentylbiphenyl) system, three different but closely related CG models, that is, 5P, 6P, and 7P, are proposed, in which the main difference is that the rigid biphenyl group is mapped into the different numbers of CG beads. Then, by analyzing and comparing the performance of these three different CG models with respect to reproducing the specific properties under investigation, the 6P model is found superior over the other two CG models. These results illustrate that small differences in the CG mapping schemes Received: November 26, 2013 Revised: March 25, 2014 Published: April 8, 2014 4647

dx.doi.org/10.1021/jp411615f | J. Phys. Chem. B 2014, 118, 4647−4660

The Journal of Physical Chemistry B

Article

monosaccharides derived at a single temperature, pressure, and concentration by force-matching method showed a good thermodynamic transferability when used in a certain range of temperature, pressure, and concentration. Louis31 found that the density-dependent CG potential will lead to transferability problems only when a resultant potential is used in state points, where the system densities differ significantly from which this potential is derived at. Therefore, although the CG potential is developed to fit one particular system property at the expenses of other properties, one can still expect that the CG potential does maintain transferability to a certain degree. To further improve “thermodynamic transferability”, scientists have taken different routes to derive the parameters for the CG potential. Krishna et al.22 proposed a thermal rescaling method and showed that the CG potential for simple Lennard−Jones (LJ) liquid and water derived by force-matching method can be transferred across different temperatures. In cases of polymer systems14,32,33 where a greater temperature range is desired, the parametrization is undertaken, separately, at different temperatures and the resulting CG potentials are fitted to include temperature dependence. By contrast, “molecular-fragment transferability” refers to the ability of a CG potential that is originally parametrized for a given molecule to perform well in other similar molecules which contain the same building blocks as these molecules. Recent work by Wang et al.24 indicated that the EF-CG (effective coarse-graining) force field developed for ionic liquid with a short side chain can be transferred to the ionic liquid with longer side chain well. In principle, the procedure to split the target molecule into small fragments (building blocks) and derive the CG nonbonded interaction parameters from isotropic liquids of these fragmental molecules is expected to maintain the transferability of CG parameters to some extent. Especially, for the complex molecules with many different CG beads, such as liquid crystal1,16 and peptides,30 such fragment-based method does not just make it simpler to derive a CG force field but also opens up the hope to reuse certain CG parameters for reoccurring building blocks. Nevertheless, the transferability of these fragment derived CG potentials to the target molecules has to be tested with great care. For instance, in our previous CG simulations,1 the CG nonbonded parameters for CG beads in the 5CB molecule were derived from isotropic liquids of fragmental molecules, that is, dodecane for C3 and C2 beads, benzene for C6 and P3 beads, and butyronitrile for CN bead. The transferability test shows that the fragment-based CG parameters for C3−C3, C2−C2, and CN−CN pairs have a good transferability as they not only can well reproduce the density and RDFs of test systems, such as octadecane and benzonitrile, respectively, but also can be extrapolated to the target 5CB molecule. Whereas the C6−C6 and P3−P3 nonbonded parameters based on liquid benzene are neither applicable to reproduce the density and structure of test systems of liquid biphenyl and hexylbiphenyl nor transferable to the rigid part (i.e., biphenyl) of LC 5CB molecule. Concurrently, there is considerable interest in optimizing the coarse-grained potential to reproduce multiproperties (e.g., several key structural and thermodynamic properties) with the aim to design a coarse-grained model which is thermodynamically and structurally consistent with an underlying atomistic simulation model.29 This consistency is of great importance for bridging between the different scales and allowing interchanges between different levels of resolution in a multiscale simulation. Recently, Rossi et al.34 reoptimized the Martini CG force field by reproducing an additional target property, that is, the

have significant consequences for the results of the CG simulations. Meanwhile, as the CG force field parametrization is a crucial step, which determines the predicting power of CG models, devising a careful parametrization guideline is thus important for reaching reliable interaction parameters. So far, several different coarse-graining techniques have been proposed to develop the force field parameters in various systems such as biomolecule2−5 and polymer6−14 systems. According to the reproduced properties, these coarse-graining methods can be classified into structure-based6−8,11−18 (which aims at reproducing a set of radial distribution functions (RDFs)), forcebased2,19−25 (monomer interactions), and thermodynamicquantities-based3−5,26−28 (density, free energy, or interface tension, respectively) ones. Furthermore, two kinds of potential forms, that is, analytical and tabulated potentials, are often chosen for the interactions between CG particles. Generally, no matter what our choice is, we should hold the same aim at achieving a simpler description of the interactions while minimizing the losses in the applicability or capability of the CG models to predict the property or properties of interest.14 In fact, fitting with high accuracy to one specific property is likely to result in a less accurate reproduction of other properties. Some recent works have indicated that the potential derived by structural-based CG algorithm may not correctly capture thermodynamic quantities, and vice versa, and the potential fitted to system thermodynamic quantities could not generate the right structural properties. Although the problem that the CG potentials derived by these coarse graining methods can correctly reproduce the target property or properties but cannot quantitatively reproduce all the properties of the system under study is now well recognized in the literature, it is not clear whether the resulting CG potential maintains transferability to some extent.29 Therefore, it is worth it to explore how far the CG potentials derived from these methods can successfully reproduce the same simulation observables under different thermodynamic conditions or for other molecules with similar chemical structures. To our best knowledge, there are still few papers studying transferability problems. In general, the term “transferability” is defined as the ability of a CG potential developed for a specific molecule under a given thermodynamic condition to predict the same observables for the same molecular system but under different thermodynamic conditions or for other molecular systems without further optimization. It can be classified into two categories, “thermodynamic transferability” and “molecular-fragment transferability” according to Engin et al.30 “Thermodynamic transferability” refers to the ability of a CG potential to reproduce the same physical quantities under the thermodynamic state points (such as temperature, pressure, density, and composition) differing from which the CG potential is explicitly parametrized for. In fact, the coarsegraining techniques mentioned above are always based on fitting one particular property, for example, density or RDF, under a given state point (a specific thermodynamic condition e.g., 300 K, 1 atm); the derived CG potentials are subsequently state-point dependent to some extent. In principle, the resulting CG potential derived in one state point cannot be used in other state points. However, the thermodynamic transferability is not as limiting as it first appears. If a resulting potential is used for state points close to the one for which it was parametrized, the transferability problems are not deemed to be important. For example, Liu et al.23 reported that the CG potential for 4648



Article

Figure 1. Coarse-grained models for (a) 6CB and (b) 8CB molecules. Circles in the atomistic 6CB and 8CB representations in the upper panels correspond to CG beads in the coarse-grained representations in the lower panels.

thermodynamic consistencies are particularly important for studying phase transitions in soft-matter systems, which are often accompanied by the variations of density and structural properties. As such, the utility of the resulting 5CB CG potentials to quantitatively modeling of nCB series is expected. In this work, we propose to study the transferability of the CG force field of the 6P model, as mentioned above, which by construction captures the structure and thermodynamic data at a given state point and can be used to investigate the phase behavior of 5CB molecules, by extending the 6P model and the associated FF to 6CB and 8CB liquid crystal molecular systems. It is worth it to note that the family of nCB shows a strong dependence of LC phase stability on the chemical details of their forming molecules, that is, nematic phase for n ≥ 5 and smectic phase for n ≥ 8. In this context, 6CB and 8CB studied herein constitute attractive test benches for the transferability problems. The remainder of this paper is organized as follows. In section two, first the CG models and interaction potentials are constructed for 6CB and 8CB based on the 6P model. Then the computational details as well as the important quantities to identify the mesophases are described. In section three, the transferability of the CG force field of the 6P model is estimated and discussed in detail by checking the nematic to isotropic (TNI) transition and nematic to smectic transition (TSN) temperatures, the molecular arrangement, as well as the diffusion coefficients in isotropic, nematic, or/and smectic phases from CG MD simulations on 6CB and 8CB systems. Finally, the prime results of this paper are summarized in the last section.

structural property of the radius of gyration, instead of just reproducing thermodynamic properties in the original Martini CG force field. Interestingly, the new CG potential had a reasonable “thermodynamic transferability” in a wide range of temperature from 350 to 600 K for polystyrene molecule. On the other hand, Allen and Rutledge35,36 derived a densitydependent implicit solvent CG potential by reproducing the RDF and solute excess chemical potential, and demonstrated this potential utilized both structural and thermodynamic properties matching and can replicate both the system RDF and excess chemical potential well when applied in mixtures at low solute concentrations.31 Further, Engin et al.30 employed the fragment-based method and developed the CG force field such that the results from CG simulations correctly reproduce both structural and thermodynamic properties of phenylalanine dipeptides. Subsequent tests showed that the CG potentials derived by the above procedure that combines both structurebased and thermodynamic quantities-based methods with fragment-based method exhibit a good “molecular-fragment transferability” to some extent. For example, when the CG potentials were used for valine−phenylalanine and isoleucine− phenylalanine systems, both the intra- and interpeptide conformational distributions as well as potential of mean force (PMF) curves and association constants are well reproduced. Similarly, in our previous work, to construct the optimal CG force field for 5CB liquid crystal molecule, we chose the aforementioned structure-based and thermodynamic quantities-based methods and used density and structural properties of the 5CB nematic phase at 300 K and 1 atm as targets for the parametrization of CG interaction potentials, and at the same time attempted to use several fragment molecular systems to derive the CG nonbonded interaction parameters. More interestingly, our results showed that the resultant CG potentials are not only thermodynamically and structurally consistent with the underlying atomistic model but also maintain “thermodynamic transferability” to some extent. For example, there is excellent consistency in the transition temperature, the transition order, and the structure and thermodynamic properties of the phases formed at different state points among these CG and atomistic molecular dynamics (MD) studies. All these findings suggest that the multiproperty parametrization route together with fragment-based method may maximize the possibility to simultaneously reproduce multiple physical properties of a given molecule or related molecules with similar chemical structures at other state points. Of course, the resulting features such as structural and

2. METHODOLOGY 2.1. Mapping Scheme. Based on the mapping scheme chosen for constructing 6P model of 5CB molecules,1 the coarse-grained models were created for both 6CB and 8CB molecules. As illustrated in Figure 1, the definition of the CG sites is given as follows. In line with 6P, the cyano head is mapped into a single CG bead of type CN and the biphenyl segment is mapped into three CG beads of P3, C6, and P3. Likewise, for the soft aliphatic tail, a mapping of three and the remaining two methylene units, respectively, into one C3 bead and one C2 bead is chosen, which is also similar to the current mapping of n-alkanes.37 All the CG beads are located at the group centers of mass in the mapping process. 2.2. Potentials. Conventionally, the total potential energy UCG is separated into bonded and nonbonded contributions. As has been stressed in the coarse-grained simulations of 5CB, the 4649



Article

Table 1. Bonded and Nonbonded Interaction Parameters for CG Simulations parameters for bond stretching interactions bond

l0 (nm)

CN−C6 CN−P3 C6−C6 P3−C6

0.335 0.242 0.424 0.308 angle CN−P3-C6 P3−C6−P3 C6−P3−C3 P3−C3−C3 C3−C3−C2

Kb (kJ·mol−1·nm−2)

fix P3−C3 fix C3−C3 fix C3−C2 fix parameters for bond angle interactions

l0 (nm)

Kb (kJ·mol−1·nm−2)

0.29 0.3392 0.29

20 000 5892 9600

Ka (kJ·mol−1·rad−2)

θ0 (deg)

360 380 95 35 30 parameters for nonbonded interactions

180 180 163 146 152

type

σ (nm)

ε (kJ·mol−1)

CN C6 P3 C3 C2

0.36 0.545 0.46 0.47 0.43

1.8963 3.5556 1.4815 1.8963 1.0667

2.3. Simulation Details. All CG MD runs were performed with the program GROMACS 4.0539 in the NTP ensemble at 1 atm and several different temperatures. The Nose-Hoover thermostat40 was employed to keep the temperature constant with a coupling time of 0.1 ps, while the Parrinello−Rahman barostat41 was used in three dimensions to control the pressure with a coupling time of 1.0 ps. Nonbonded interactions were truncated by the twin-range cutoff method with a short cutoff of 1.2 nm and a long one of 1.5 nm. During all runs for 6CB and 8CB, the rigid bond lengths were kept fixed at their equilibrium value with the Lincs algorithm and a time step of 10 fs was employed. In this work, all the simulation were performed on a system including N = 2048 6CB or 8CB molecules with ordered starting configurations. Such configurations were obtained as follows: First, all molecules were placed with their long axes pointing along the x axis of the simulation box; then the system was heated up to 330 K and equilibrated until it evolved to a disordered state with the orientational order parameter P2 less than 0.1; finally, the system was cooled down to 300 K for about 200 ns to get an ordered starting configuration with P2 = 0.5 for 6CB and P2 = 0.55 for 8CB for the later runs at different temperatures. For each temperature, the equilibration of these runs was accessed by monitoring total energies, densities (ρ), and orientational order parameters. Because of the long time scale involved, particularly near the TNI or TSN, the simulations at each temperature were run sufficiently long to reach equilibration, and then additional production runs of at least 200 ns were used in order to yield reliable results. To monitor the evolution of the system, the internal energy, the bulk density, and the orientational order parameter P2 have been calculated. P2 is defined as the largest eigenvalue of the Saupe ordering matrix Q42

intramolecular bonded interactions are derived such that the local conformational distributions of the molecules are reproduced in the CG MD simulations. Keeping in line with the bonded interaction potentials of the 6P model, only bond stretching and bond bending interactions are considered, and two harmonic potentials are used for them respectively in each individual CG 6CB and CG 8CB molecules., Vbond = 0.5Kb(l − l0)2

(1)

Vangle = 0.5K a(θ − θ0)2

(2)

where Ka and Kb are the CG force constants, respectively, and l0 and θ0 are the CG equilibrium bond length and bond angle, respectively. In the present work, the FF parameters for CG bond stretching and bond bending interactions are directly taken from 6P model except for the CG bond stretching interaction of C3−C3 and bond bending interactions of P3− C3−C3 and C3−C3−C2. These three interaction parameters have been derived from sampling the corresponding distributions of the atomistic simulation of the same molecular system using the approach introduced in section 2.2.1 of ref 1. As for the CG nonbonded interaction potential, following the 6P model, a common functional form of the soft LJ 9-6 potential is selected for nonbonded interactions in the CG 6CB or CG 8CB simulations U9‐6(r ) =

bond

9 ⎛ σ ⎞6 ⎞ 27 ⎛⎜⎛ σ ⎟⎞ ε⎜ −⎜ ⎟⎟ ⎝r⎠ ⎠ 4 ⎝⎝ r ⎠

(3)

where the parameters ε and σ are the minimum energy and the distance at U = 0, respectively. The standard Lorentz− Berthelot38 mixing rules of εij = (εiεj)1/2 and σij = (σi + σj)/2 were used for calculating the LJ parameters between different types of CG beads. The nonbonded parameters we adopted here are also taken from 6P model, which were originally developed to reproduce multiproperties including density and radial distribution functions as well as nematic phase of 5CB systems at 300 K and 1 atm. All the bonded and nonbonded interaction parameters are listed in Table 1.

Q ab =

3 1 uaub − δab 2 2

(4)

where the mean value ⟨...⟩ is obtained by averaging on all the molecules and ua (a = x, y, z) is the vector of molecular long axis defined on the basis of the CN−C6 vector (Figure 1). The 4650



Article

3. RESULTS AND DISCUSSION As the main concern of the present paper is to study the transferability of our previously developed FF parameters for the coarse-grained 6P model of 5CB molecules, it would be highly desirable to study whether 6P model parameters can quantitatively reproduce the phase behavior (i.e., TNI and TSN), some structural properties, as well as diffusion coefficients of the 6CB and 8CB liquid crystal systems. 3.1. 6CB Molecular System. First, we test whether the CG 6CB model derived by directly extending the 6P model of 5CB can reproduce the nematic phase. We performed a series of CG MD simulations on a 6CB system at temperatures between 295 and 310 K from an identical orientationally ordered starting configuration with the orientational order parameter P2 = 0.5, and the resulting P2 values are shown in Figure 2a. At the

director n is identified as the eigenvector corresponding to P2. The value of P2 will be close to 0 in the isotropic phase, but tends to 1 in a highly orientational ordered phase. In our simulation, we could distinguish the isotropic phase from the nematic phase by looking at this order parameter. Additionally, because of the finite system size, the value of P2 in the isotropic phase is small but nonzero. In order to fully characterize the phase, the detailed structures of the LC phases were also investigated by the various distribution functions. One quantity is the orientational pair correlation function, which measures the decay of intermolecular angular correlations as a function of their distance and can thus be used to quantify the orientational order of a phase. Two orientational correlation functions are computed as g1(r ) = ⟨cos θij(r )⟩

(5)

3 1 (cos θij(r ))2 − 2 2

g 2 (r ) =

(6)

where θij(r) is the angle between the long axes (as defined above) of two molecules i and j, and r is the distance between the reference sites on two different molecules, that is, the center of mass of cyano group. It is worth mentioning that g1(r) can give information on the local polar ordering, while g2(r) becomes long-ranged in the NI transition and its ultimate value at large r tends to ⟨P2⟩2. The RDF from atomistic simulations is used as one of target quantities to derive the CG interaction potentials. We calculate the RDF between a site k in molecule i and a site m in molecule j as g (r ) =

1 ⟨δ(|rik − rmj|−r )⟩ 4πr 2ρN

Figure 2. (a) Time evolution of orientational order parameter for the coarse-grained simulations of 6CB systems at several temperatures. (b) Averaged order parameter P2 and density ρ as a function of temperature for CG 6CB systems and the underlying atomistic systems, together with several sets of the experimental densities. The blue triangles represent experimental values,43 and the vertical dash line indicates the experimental TNI..

temperatures of 295 and 300 K, the systems remain orientationally ordered, since the orientational order parameters over the last 200 ns are about 0.55 at T = 295 K and 0.48 at T = 300 K, which are within the typical P2 range from 0.4 to 0.6 for the nematic phase.44 Particularly, the simulated P2 value of 0.55 at 295 K is very close to the experimental value of 0.51.45 However, as the temperature is increased to 305 and 310 K, P2 goes down to a low value, that is, ⟨P2⟩ over the last 200 ns is about 0.2 for T = 305 K and 0.05 for T = 310 K. Obviously, these two values are so small that 6CB systems at 305 and 310 K can be considered as in the isotropic phase, since the order parameter of an isotropic phase is usually below 0.4.44 Hence, we can tentatively locate the NI transition temperature of CG 6CB systems between 300 and 305 K, which is very close to the experimental value of 302 K,46 indicating the force field parameters originally developed for CG 5CB model can reproduce the NI transition temperature of 6CB molecular system very well. Meanwhile, the presence of the NI phase transition can also be revealed by the temperature-dependent averaged order parameter and/or density. As presented in Figure 2b, the CG simulated ⟨P2⟩ exhibits a jump between the temperature of 300 and 305 K, in well agreement with atomistic simulation results, whereas the density profiles from both CG and atomistic simulations do not clearly show such a sharp discontinuity around the NI transition. Similar findings have been reported by us from CG 5CB simulations as well as atomistic 5CB MD studies, and all are in agreement with the weak first-order nature of the NI transition;47 that is, NI transition sometimes cannot be absolutely discrete by the densities. More importantly, the resulting densities at 295 and

(7)

where the mean value ⟨...⟩ is a double average over all molecules and all time origins, and ρN is the number density. The onset of the smectic order is monitored by the onedimensional positional order parameter τ, which is defined as τ=

1 N

N

⎛

∑ exp⎜⎝2πi k=1

rk ·n ⎞ ⎟ d ⎠

(8)

where N is the total number of molecules, rk is the coordinate of the center of mass of the kth molecule, and d is the smectic layer spacing. For the smectic-A phase, the layer normal coincides with the director of the system. As the layer spacing is not known a priori, we optimized in order to give the maximum value of τ. τ vanishes in the isotropic and nematic phases, is nonzero in the smectic phase, and is unity in the limit of complete one-dimensional translational order. To further explore the reliability of our FF, we also investigate one important aspect of the dynamical behavior, namely, the translation diffusion. The diffusion coefficients Dα are computed from ⎛1⎞ Dα = lim ⎜ ⎟ t →∞⎝ 2t ⎠

N

∑ [R jα(t ) − R jα(0)]2 j=1

(9)

Rαj

where α = x, y, z. is the position of the center of mass of molecule j projected onto the reference axis. Note that to calculate Dα additional 50 ns simulations under NVT ensemble were performed just after the long simulations under NPT. 4651



Article

in line with our previous results on 5CB molecules.1 In particular, the discrepancies in g1(r) peak width and height between atomistic and coarse-grained models become stronger when the simulation temperature is further from which this CG FF is derived at. For example, as typically illustrated in Figure 3a, the g1(r) peak heights at 300 K are 0.83 from CG simulations and 0.88 from atomistic simulations, while they turn to be 0.74 and 0.84 at 305 K. As coarse-graining the groups of atoms into spherical CG beads with radically symmetric CG potentials often bears the problem that the link to the underlying structure and the reproduction of real local packing effect may become weak, certain discrepancies at the level of local polar ordering between the atomistic and the CG description are not unsurprising. Specifically, instead of tuning the CG FF parameters to reproduce the orientational correlation function, our CG force field of 6P model is obtained by reproducing the target bulk density and RDF; thus, it seems difficult to precisely reproduce antiparallel arrangement details as that are present in the original atomistic system. Meanwhile, electrostatic interactions are not explicitly modeled, which may also limit the ability of our CG force field of the 6P model in an accurate reproduction of the short-range dipole correlation as the underlying atomistic model. Moreover, the enhanced discrepancy with more deviation of the simulation temperature from the one for which the CG FF is explicitly parametrized still indicates that the CG force field of the 6P model may have a limited range of transferability in quantitative reproduction of some local structural properties at other state points. Consequently, caution must be used if we extend the straightforward adoption of the CG FF parameters over other state points for which they are not explicitly parametrized, especially when quantitatively studying the local structural properties. Additionally, we should mention that this degree of discrepancy in g1(r) peak width and height, however, does not disrupt the formation of antiparallel arrangements between the nearest neighboring 6CB molecules. Also, as can be seen that in Figure 3a, the positions of the negative peaks in g1(r) curves from the CG simulations are almost coincident to ones from the atomistic simulations. Therefore, regardless of some discrepancies in g1(r) peak width and height, the CG force field parameters of 5CB molecules at 300 K maintain transferability to some extent and are able to predict the antiparallel arrangements of 6CB molecules at short distance in a temperature range we studied. At the same time, the RDFs for the pairs of other CG sites of C3−C3, C6−C6, and P3−P3 obtained from both atomistic and CG simulations on the 6CB at T = 300 and 305 K are calculated and given in Figure S1 in the Supporting Information. It should be noted that the positions of the RDF peaks, especially of the first one, from CG simulations for all types of same CG bead pairs except for C6−C6 closely match with those from atomistic simulations. In particular, although the electrostatic interactions which are important for atomistic simulations are ignored in our CG simulations, the position of the first peak of RDFs for the CN−CN pair still matches to each other. A relatively large deviation occurs for the C6 bead, whose first peaks are located at 0.6 and 0.57 nm for CG and atomistic simulations, respectively. As already mentioned before, this discrepancy at the level of local packing between the atomistic and the CG description is not unexpected. Thus, the observation that resulting positions of the first RDF peaks in our CG simulations are comparable with those of atomistic simulations further confirms that our CG

310 K are 1002 and 985 g/mm3, very close to the experimental values of 1020 and 995 g/mm3 as well as the atomistic simulated results of 1037 and 1025 g/mm3, respectively, wherein the corresponding density underestimations are less than 2% when compared with experimental densities and about 3% in comparison with atomistic simulation estimates. This illustrates again that although the current CG force field is parametrized for 5CB molecules at a temperature of 300 K, it maintains transferability to some extent and is valid for predicting the densities of 6CB systems in a temperature range we studied. After the consistency in TNI and the IN transition order among these CG and atomistic MD studies on 6CB has been verified, it is of particularly interest to study the arrangement of 6CB molecules in both nematic and isotropic phases by the orientational correlation function g1(r) and the radial distribution functions g(r), since as a homologues of 5CB 6CB molecules also prefer antiparallel arrangements at short distance due to the electrostatic interactions. As seen from Figure 3a, the g1(r) data from the CG simulations are first

Figure 3. (a) Orientational correlation function g1(r) and (b) radial distribution function as a function of distance between the centers of mass of the CN groups from atomistic (labeled as AT) MD and between CN beads from CG (labeled as CG) MD simulations at T = 300 and 305 K.

negative up to 0.5 nm at T = 300 and 305 K, which is the quite similar to the atomistic simulation results. Moreover, for the CN−CN bead pair, the position of the first peak g(r) from the CG simulations is about 0.39 nm at T = 300 and 305 K, indicating that the most favorable geometrical arrangement between the nearest neighboring molecules in nematic (300 K) and isotropic (305 K) phases is also antiparallel in the CG simulations. Similar antiparallel arrangements have been found in CG simulation studies of 5CB molecule system.1 As addressed in our previous work,1 this antiparallel arrangement is a result of the coulomb interactions, which are essentially determined by the presence of the polar CN groups in 6CB molecules. Although in our CG 6CB model the dipolar interactions are not considered explicitly, the adopted CG force field was originally constructed in a way to reproduce the nematic structure features of the underlying atomistic model for 5CB at 300 K, thus the local structure information associated with the electrostatic interactions is already incorporated into the nonbonded interaction parameters of CN beads to some extent. Obviously, the occurrence of antiparallel arrangements demonstrates the validity of CG FF parameters by straightforward adoption of the CG 6P model of 5CB molecules in accounting for the short-rang dipole correlation. However, it is worth noting that not only the peak width of g1(r) from CG simulations is thinner than the corresponding atomistic one but its negative peak value from CG simulations is also smaller than the one from atomistic simulations, and this observation is also 4652



Article

force field of 6P model can maintain the true excluded volume “shape” of CG beads in 6CB molecules as close as possible and is applicable to properly reproduce the local packing effect to some extent. On the other hand, we found from the RDFs in Figure 3b as well as in Figure S1 in the Supporting Information that most RDFs in our CG simulations display a relatively higher first peak than those obtained from the corresponding atomistic simulations. Although, due to the necessary simplifications inherent in the creation process of a CG model, we are unable to completely reproduce the correct packing with the CG model; the stronger peak heights in the CG model still reflect an effect due to the use of a softer LJ 9-6 nonbonded interaction potential, which enhances the overlapping of CG beads and thereby affects the local packing.1 In general, the soft interaction potential could lower the effective repulsion between CG beads and reduce the effective density of the system. Therefore, it is conceivable that this high packing in the CG model is in line with the preceding CG results, wherein the density is found to be 2% and 3% below the experimental and atomistic simulated data, respectively. To yield a better representation of the real local packing with a larger CG particle, a nonbonded potential with the steeper repulsion term is required.1 However, it should be noted that this degree of discrepancy in local packing is not a critical issue for the 6CB liquid crystal system. In addition to the thermodynamic and structural properties, we also extend the validation procedure on the transferability of the 5CB CG FF to the dynamic properties of the 6CB. It is well-known that the anisotropic diffusion of the nematic phase has been the subject of many experimental as well as theoretical investigations. So it is quite interesting to study this special property in a CG 6CB model system. Then the diffusion coefficient for the 6CB nematic phase at 295 K is calculated. As shown in Table 2, D∥ (parallel to the director) is larger than D⊥

Information when compared to the underlying atomistic representation, slows down the translational movements. Moreover, the levels of coarse-graining in our 6CB and 8CB systems wherein only two or three atoms are mapped into one CG bead are relatively low if compared to the other systems, such as polymer and biomolecules wherein 10 or more atoms are mapped into one CG bead. This also could decrease the speedup in the dynamics of CG simulation. Overall, our CG 6CB model can reproduce not only the NI phase transition temperature, the experimental density, the neighboring molecular antiparallel arrangement, and the RDF positions, but also the diffusion properties of 6CB nematic phase. Additionally, it is worth it to stress that since only one extra methylene is added to the alkyl chain of 6CB molecule when compared with the 5CB molecule, the actual difference in TNI between 5CB and 6CB is extremely small, of about 5 K. Nevertheless, such a subtle deviation is also distinguished from our simulations on CG 6CB and 5CB models; that is, the resulting TNI is 305−310 K for 5CB 6P model systems and is between 300 and 305 K for 6CB CG model systems. Apparently, all these not just validate the rationality of our constructed CG 6CB model, but more importantly also confirm the transferability of the CG force fields developed for 5CB molecules to some extent. Moreover, it deserves noting again that to construct this CG force field we combine the thermodynamic quantities-based and structure-based methods and devise a careful guideline to parametrize the CG interaction potentials in a systematic way, and at same time attempt to use several fragment molecular systems to derive the CG nonbonded interaction parameters. Precisely for these reasons, the derived CG potentials could both properly maintain the 5CB CG model thermodynamically and structurally consistent with the underlying atomistic model1 and are transferable to 6CB liquid crystal systems across different thermodynamic conditions we studied. 3.2. 8CB Molecular System. We should stress again that different from 5CB and 6CB systems which only have nematic mesophase, two mesophases, that is, nematic and smectic phases, have been observed experimentally in 8CB.54,55 In view of such fairly rich phase behavior as well as important applications of 8CB, 8CB presents a much more interesting test system to access the above addressed transferability issue. A good transferability will be verified if CG MD simulations with our developed 5CB CG FF can reproduce the phase behavior, some key structural properties, and diffusion behavior of both nematic and smectic phases of 8CB. 3.2.1. Phase Transition. First, we check whether the CG 8CB model derived by straightforward adoption of 6P model parameters can reproduce stable nematic and stable smectic phases near the experimental range of existence. A series of CG MD runs have been performed at several different temperatures from an identical initial configuration with an orientational order parameter of P2 = 0.55 and positional order parameter τ < 0.1 (in other words, the starting configuration is in a nematic state). The resulting P2 values are shown in Figure 4a, the orientational correlation function g2(r) in Figure 4b, and typical snapshots in Figure 5. At 330 and 320 K, the 8CB system loses its orientational order, as P2 decreases remarkably and ⟨P2⟩ over the last 200 ns is less than 0.1. Meanwhile, g2(r) decreases to near zero at 320 K and above, showing the feature of isotropic fluid. These, together with the random orientational arrangements of 8CB molecules in the snapshot of Figure 5a, indicate that, as the

Table 2. Comparison of Diffusion Coefficients between CG Simulations and Experiment for the 6CB Molecular System translational diffusion coefficients

CG simulation (295 K)

experiment48 (296 K)

D∥ (10−11 m2/s) D⊥ (10−11 m2/s)

4.1 2.1

4.5 none

(perpendicular to the director), yielding the expected diffusion anisotropy with the ratio D∥/D⊥= 1.95. However, in contrast to the mostly observed dynamics speedup in CG simulations, the resultant value of D∥ is about 4.1 × 10−11 m2/s at 295 K, which is close to the available experimental value48 of 4.5 × 10−11m2/s at a nearby temperature, that is, 296 K. Note that Izvekov and Voth49 have addressed recently that to achieve a more faithful representation of the real dynamics in the CG simulation, the friction and noise terms are needed to be introduced to compensate the missing degrees of freedom during the coarsegraining mapping. Along this line, the dynamics in the coarsegrained simulations is found to be corrected by solving the generalized Langevin equation50 or by using the Langevin thermostat51 or DPD thermostat.52,53 While, without such a correction in our CG simulation, the nonaccelerated dynamics could be caused by several reasons. The most important one we think lies in the excluded volume effect, wherein the presence of a big spherical C6 bead in the 6P model, as evidenced by a higher nearest-neighbor peak at a relatively larger position in the RDF of Figures S1d and S2c from the Supporting 4653



Article

the temperature is further decreased to 290 K and below, P2 increases and ⟨P2⟩ over the last 200 ns are all larger than 0.5, that is, about 0.62 for 290 K and 0.78 for 280 K. Particularly, at 280 K, as shown in the snapshot in Figure 5c, a layer-like structure appears. From the data we gathered at these temperatures, we come to the conclusion that at T ≤ 315 K 8CB systems prefer to maintain in the orientationally ordered states, which is further supported by the asymptotic behavior of g2(r) to a stable value of ⟨P2⟩2 at large distance at 280, 290, and 315 K as shown in Figure 4b. To further quantitatively judge whether the resulting orientationally ordered phase is a smectic or a nematic phase, the translational order parameter τ is used to characterize the evidence of smectic phase. As shown in Figure 4c, at 280 K, τ has an apparent broad peak with the maximum value of 0.45 at the distance about 2.0 nm, suggesting again that the orientationally ordered phase at 280 K possesses some layer-like structures. While at 290, 315, and 320 K, there are no peaks with heights above 0.1, indicating that no positional order is present. Hence, the 5CB CG FF reproduces not only the nematic phase but also the smectic phase of 8CB. From these results, we can locate the NI transition temperature from our CG 8CB model between 315 and 320 K, which is close to our atomistic simulation result of TNI = 305−310 K,56 and both simulated data are in good agreement with the experimental value of 313.5 K.55 The smectic−nematic transition temperature (TSN) from our CG 8CB model is between 280 and 290 K, while our atomistic simulation result of TSN = 300−305 K56 and the corresponding experimental measure is 306.5 K.54 Meanwhile, the presence of the NI and NS phase transitions can be clearly revealed by jumps in the temperature-dependent ⟨P2⟩ in Figure 4d. Moreover, similar to the results from previous atomistic simulation of 8CB56 as well as 6CB CG simulated results in Figure 2b, there are no sharp discontinuities in the density curves, indicating a weak first-order nature of the NI transition,47 and also in accordance with the expectation that the N−S transition should become more weakly first order for elongated molecules. Interestingly, the density of our CG simulated result for the 8CB nematic phase at 310 K is slightly lower than the available experimental measurement and atomistic simulated result; that is, the underestimations are less than 1% and 2%, respectively, indicating that the density in the 8CB nematic phase can be reproduced well to some extent. Meanwhile, it can be seen that the resulting density of CG simulation at 300 K is also comparable with the corresponding experimental value and atomistic simulation data with an underestimation about 1% and 2%, respectively, but CG MD simulations cannot exactly reproduce the smectic phase as observed in atomistic simulations and experiments at 300 K and 1 atm and thus lead to the incorrect behavior. In short, although the current CG force field is parametrized for the nematic 5CB phase at 300 K and 1 atm, it can reproduce TNI up to a satisfactory level of accuracy for both 6CB and 8CB. However, it does not perform sufficiently well for predicting the smectic−nematic transition temperature; that is, the difference between CG simulated and experimental values of TSN is about 20 K, beyond the temperature scan step of 10 K in CG simulations. These indicate that, in spite of this CG FF incorporating both thermodynamic and local structural information on the nematic 5CB phase at 300 K, it is to some extent transferable to the 6CB and 8CB nematic phases. While this CG FF is not well suited to reproduce the smectic

Figure 4. (a) P2 as a function of time for several different temperatures. (b) g2(r) as a function of distance between CN beads at four different temperatures of 280, 290, 315, and 320 K. (c) The translational order τ averaged over the last 100 ns of CG simulation is calculated as a function of distance between CN beads at four different temperatures of 280, 290, 300, and 330 K. (d) ⟨P2⟩ and densities from CG simulations (black line) as a function of temperature. For comparison, densities from atomistic simulations56 (red line) and the available experimental density values57 (blue triangles) are plotted. The two vertical dashed lines (black) denote the experimental TSN (306.5 K) and TNI (313.5 K). The other two vertical dotted lines (red) represent the TSN (280−290 K) and TNI (315−320 K) from CG simulation.

Figure 5. Snapshots of last configurations from CG 8CB simulations at (a) 330 K, (b) 310 K, and (c) (280 K). In all these pictures, the CG 8CB molecule is presented by the rigid part which contains CN (blue), P3 (yellow), C6 (red), and P3 (yellow) beads interconnected by bonds.

temperature increases to 320 K and above, the stable state for the 8CB system is an isotropic phase. In contrast, at the temperature T = 315 K, P2 first drops slowly and then grows up at about 200 ns and thereafter fluctuates between 0.4 and 0.55 for more than 250 ns. Similarly, at 310 K, the 8CB systems remain orientationally ordered, as ⟨P2⟩ over the last 200 ns is about 0.48 and 8CB molecules lie parallel to a common direction as clearly seen from the snapshot in Figure 5b. When 4654



Article

interactions and enhancement of the anisotropy of molecular shape by increasing the density or lowering temperature work in concert to postpone the smectic phase to lower temperatures at which 8CB systems have higher densities and higher effective length-to-width ratios. Lastly, we should keep in mind that electrostatic interaction is not explicitly modeled, which may also limit the scope of our CG force field of 6P model to precisely reproduce the experimental TNS as the underlying atomistic model. 3.2.2. Molecular Arrangement. After the liquid crystal phases have been determined, it is interesting to study the molecular arrangement. The orientational correlation function g1(r) and the radial distribution function g(r) are calculated as a function of distance between CN beads in order to characterize the local arrangement of CG 8CB molecules in two liquid crystal phases, i.e., smectic (280 K) and nematic (305 K). As presented in Figure 6a, all g1(r) curves from CG and atomistic

phase of the 8CB system at 300k and 1 atm and the smectic phase is postponed to lower temperatures at which 8CB systems have higher densities and higher effective length-towidth ratios. Indeed, the lack of transferability cross the smectic phase is not a desirable situation. In general, as one of the inevitable consequences of the coarse-graining technique, the necessary simplifications inherent in the CG process make it impossible to quantitatively reproduce all the properties of a real system with a coarse-grain model. Therefore, it seems possible that a CG model, derived on the basis of a given scheme and some approximations, does not display the N−S phase transition at the very same temperature and pressure as the underlying atomistic model.58 Meanwhile, as mentioned before, by construction this 5CB CG FF inevitably incorporates the detailed structural information regarding the specific nematic phase of 5CB at 300 K and 1 atm, which may not allow for a locally very tight packing of 8CB molecules and thus fails to reproduce the smectic phase at 300 K and 1 atm. Increasing density can enhance the local packing of the molecules. From this, we can infer that TSN for the CG model will be at lower temperatures than for the chemically realistic one. Besides these possible explanations, we can also tentatively interpret the limitation of current 5CB CG FF in reproducing the semctic phase points in light of the CG interaction potentials. As addressed before, the soft nonbonded potential used here lowers the effective repulsion between CG beads and reduces the effective density of the system. It has long been recognized that excluded volume effect has a large effect on the formation of LC mesophases. Thus, it is not unexpected that small density underestimation of CG simulated 8CB at 300 K and 1 atm may even cause the smectic phase not to be reproduced at the CG level. In addition to excluded volume, anisotropy in the shape of a molecule is a fundamental ingredient necessary for the formation of liquid crystal phases according to Onsager’s theory.59 Generically, long rods order themselves at lower densities than short rods. It should be noted that the actual length-to-width ratio of mesogens is a decreasing function of flexibility and temperature.60 Decreasing the stiffness of the mesogens reduces the anisotropy of molecular shape and thereby shifts the liquid crystal phases to the lower temperature or the higher density end of the phase diagram.60,61 In the current 5CB CG FF, we excluded all intrachain bonded interactions beyond the distance of three CG bonds (i.e., torsions between any quadruple of beads), which may affect the overall stiffness of the chain though this effect may be relatively small because the lincs method was used to constraint bond lengths in the rigid biphenyl unit. In other words, despite the 8CB molecules being able to retain their basic rodlike structure at a CG level, a removal of intramolecular barriers against rotation will more or less increase internal disorder; the resulting some flexibility in the CG model may lead to a decrease of shape anisotropy of the biphenyl segment if compared with the corresponding atomistic measurement. This makes formation of the smectic phase difficult at 300 K and 1 atm and thus favors the nematic phase when using current 5CB CG FF. As just mentioned, one can enhance the effective length-to-width ratio by decreasing the temperature. Also, mesogens become “longer” and “thinner” as density is increased. From this, we conclude that smectic phases are therefore pushed to the high density or low temperature end of the phase diagram. Of course, the two explanations based on using soft nonbonded potential and excluding torsion barriers are correlated: strengthening excluded volume

Figure 6. (a) Orientational correlation function g1(r) and (b) radial distribution function as a function of distance of the CN groups from atomistic (labeled as AT) MD and between CN beads from CG (labeled as CG) MD simulations at T = 280 K (smectic) and 305 K (nematic).

simulations at T = 280 K (black line) and 305 K (red line) first show negative peaks at a short distance that is very close to the first peak positions of the corresponding g(r) curves. This means that, like the underlying atomistic models, the most favorable geometrical arrangement between the nearest neighboring 8CB molecules in nematic (305 K) and smectic (280 K) phases is also antiparallel in the CG simulations. Obviously, the occurrence of antiparallel arrangements demonstrates again that information associated with the short-range dipole correlation was incorporated in the 5CB CG FF. Also, like what we found in 5CB and 6CB molecules, some discrepancies in g1(r) between CG and atomistic simulations have been observed. For example, as typically illustrated in Figure 3a, the g1(r) negative peak height at T = 305 K is 0.72 from CG simulations, which is smaller than the atomistic one of 0.83. Additionally, the peak widths of g1(r) from CG simulations at T = 280 and 305 K are thinner than the corresponding atomistic ones. As addressed for 6CB, the reason for these discrepancies may come from the simplifications adopted in our coarse-graining process; for example, we chose symmetric CG potentials which do not explicitly include electrostatic interactions and are parametrized by reproducing the target bulk density and RDF. Furthermore, we note that for the smectic phase (T = 280 K) a strong positive peak shows up in the g1(r) curve at a distance of about 0.6 nm and is followed by several progressively reduced periodic peaks. In contrast, the g1(r) data for the underlying atomistic model converges directly to zero beyond a weak positive peak at a similar distance of 0.6 nm. Such a deviation indicates again that the CG force field of 6P model may have a limited transferability in quantitative 4655



Article

a single-layer-like rather than partially interdigitated bilayer structure is generated in the CG simulation at T = 280 K. According to our previous atomistic work, the polar sublayers in the 8CB Sm-Ad phase can be detected easily by plotting the distribution function of the “up” or “down” molecules as a function of the position of R along the director, wherein the “up” or “down” molecule is identified by the angle between the molecular long axes and the director either smaller or larger than 90°, respectively. To further explore this single-layer-like structure, the distribution function of the “up” or “down” molecules was also calculated and presented in Figure 7b. Interestingly, we see that the “up” and “down” peaks oscillate synchronically and nearly overlap each other, and the averaged distance between two consecutive peaks is also about 2 nm, very close to the above derived smectic layer spacing. Obviously, within each smectic layer, about half of 8CB molecules orient to one direction, while the rest of molecules orient to the opposite. The snapshot of the 8CB smectic phase from CG simulations at 280 K confirms this arrangement. As shown in Figure 5c, contrary to the atomistically simulated 8CB Sm-Ad phase in which the blue CN sites are found to mostly distribute within the interdigitated region between the two sublayers, the blue CN sites are evenly located at the upper or lower positions within each single-layer region. Based on the above analysis together with the observed antiparallel arrangement at short distance, we can conclude that 8CB in CG simulations prefers to form a single-layer-like structure and within each layer the nearest neighbor molecules are antiparallel to each other. For a clear comparison, we draw qualitative pictures of the most probable molecular arrangements for the single-layer-like smectic phase as obtained from CG simulation and the partial bilayer-like smectic phase from atomistic simulation in Figure 8a and b, respectively. It should be

reproduction of some local structural properties at other unlike state points, that is, smectic phase. In addition, we should mention that regardless of the discrepancies in g1(r) peak width and height, the CG FF originally developed for 5CB molecules at 300 K maintains the transferability to some extent and is good enough to predict the antiparallel arrangements of 8CB molecules at short distance. However, more attention should be paid if we extend the straightforward adoption of the CG FF parameters over other state points (e.g, smectic phase) for which they are not explicitly parametrized, especially when quantitatively studying the local structural properties. Also, as can be seen from the RDF for CN in Figure 3b as well as from the RDFs for the pairs of other CG sites of C2− C2, C3−C3, C6−C6, and P3−P3 in Figure S2 in the Supporting Information, we found that they exhibit the very similar behavior as 5CB and 6CB molecules. The positions of the RDF peaks, especially of the first one, from CG simulations for all types of same CG bead pairs except for C6−C6 closely match with those from atomistic simulations, while the heights of the first RDF peaks from CG simulations are mostly higher than those obtained from the corresponding atomistic simulations. As addressed before, besides the inherent simplifications and necessary scheme in the creation process of CG FF, the soft nonbonded potential used here may also cause us unable to completely reproduce the correct packing with CG model. Furthermore, our previous atomistic simulations56 showed that, in marked contrast to the standard smectic A (SmA) phases, the 8CB smectic phase is a partial bilayer smectic (SmAd) one, wherein two sublayers, with opposite polarization, partially interdigitate to form smectic layers. As no charge or electrostatic interaction term was explicitly considered in our 5CB CG force field, a more interesting test on the transferability of the 5CB CG FF is to verify the model’s capability to reproduce the 8CB Sm-Ad phase. First, the radial distribution function of the CN bead along the director direction g∥(r) was calculated at several temperatures and is shown in Figure 7a. Obviously, only at T = 280 K, g∥(r) has the

Figure 8. Most probable molecular arrangements (a) in the singlelayer-like 8CB smectic phase from CG simulation and (b) in the partial bilayer-like 8CB smectic phase from atomistic simulation.56 Each arrow schematically represents one 8CB molecule and the arrow direction points to the C6−CN vector. The vertical dashed black lines are used to define the boundary of one smectic layer, and the vertical dashed green lines mark the positions of the center of mass of the two sublayers within each smectic layer.

Figure 7. (a) Radial distribution function of CN bead along the director in CG simulations of 8CB systems at several different temperatures. (b) Projection of center of mass of molecules along the director for a CG 8CB system simulated at 280 K. The symbol “ up” or ”down” means the angle between the director vector and principal axis is smaller or larger than 90°, respectively.

noted that the single-layer-like smectic structure can be taken as a completely interdigitated instead of partially interdigitated bilayer structure. Apparently, although the used 5CB CG force field is able to reproduce the antiparallel arrangements of the nearest neighboring 8CB molecules up to a satisfactory level of accuracy, it fails to reproduce the partially interdigitated bilayerlike smectic structure. Of course, the inherent simplifications and necessary scheme in the creation process of CG FF (for example, we chose symmetric CG potentials which do not explicitly include electrostatic interactions and are parametrized by reproducing the target bulk density and RDF of the nematic

well-defined oscillatory form, showing the tendency of molecules to arrange into layers. The averaged distance between two consecutive peaks, which corresponds to the smectic layer spacing (dsm), is about 2 nm, agrees well with the value of d at which τ plot reaches its maximum value in Figure 4b. Note that dsm from our previous atomistic simulation56 at T = 300 K is found to be 2.85 and 2.8 nm at 280 K, while the molecular length of 8CB is about 1.8 nm. Thus, we expect that 4656



Article

most likely reason for this nonaccelerated dynamics lies in the excluded volume effect of the big spherical C6 bead in the 6P model, which brings larger frictional resistance to the translational movements. Besides, the relatively less level of coarse-graining in our CG model could also contribute the nonaccelerated dynamics. Interestingly, the diffusion in the smectic phase obtained from CG simulations is still nematic-like in the sense that D∥ > D⊥. An analogous behavior has recently been revealed by NMR measurement,64 Prampolini’s simulation,62 and our previous atomistic simulation.56 At the same time, the ratio of D∥/D⊥ for the CG 8CB smectic phase at 280 K is 1.75, which can be well compared with atomistic simulation result of 1.73 at 280 K. As addressed in our previous work, it is not so unexpected that the inequality D∥> D⊥ remains on entering the 8CB smectic phase, because the local molecular environment in the nematic and smectic phases is not so dissimilar, especially when the NS transition is relatively weak or continuous. In a word, all these results indicate that, in spite of this CG FF incorporating both thermodynamic and local structural information on the nematic 5CB phase at 300 K, it is transferable to study the dynamic properties of 8CB isotropic, nematic, and smectic phases.

5CB phase at 300 K and 1 atm) may make it impossible for the CG 8CB model to precisely reproduce the bilayer-like smectic structures as the underlying atomistic model. For example, as reported by Prampolini and De Gaetani,62 when removing point charge interactions from the force field during the atomistic simulation of 8CB, the partial bilayer Sm-Ad was collapsed into the single-layer smectic phase. They therefore claimed that although the electrostatic interaction energy is less than 10% in the overall interaction energies of the smectic phase, it dominates the formation of the partial bilayer structure in 8CB smectic phases. Meanwhile, this 5CB CG FF inevitably incorporates the detailed structural information regarding the nematic phase of 5CB at 300 K and 1 atm, which may not allow for reproducing the specific bilayer-like smectic structure. Additionally, as addressed before, the soft LJ 9−6 nonbonded potential used here as well as the removal of intramolecular barriers against rotation in our CG simulation will more or less increase the molecular flexibility of CG model, which can facilitate the interpenetration of smectic layers.63 3.2.3. Diffusion Coefficients. As already mentioned, comparing the CG simulated dynamic properties with the corresponding atomistic simulation results and available experimental measurements64 is also of importance for validation on the transferability of the employed 5CB CG FF. The diffusion coefficients for the CG 8CB systems in three different phases, that is, isotropic, nematic, and smectic phases, have been calculated, and the typical results along with the corresponding experimental values are listed in Table 3. Similar

4. CONCLUSION In this paper, the transferability of the coarse-grained force field originally developed for 5CB liquid crystal molecule was investigated by its homologues 6CB and 8CB molecules. Overall, not only the experimental densities, the RDF positions, and the nearest neighboring molecular arrangement, but also the unique LC mesophases as well as the nematic−isotropic phase transition temperature of 6CB and 8CB were reproduced. Apparently, although the current 5CB CG force field is parametrized for the nematic 5CB phase at 300 K and 1 atm, it exhibits a good “thermodynamic transferability” and “molecular-fragment transferability” to some extent, which further validates the rationality and effectiveness of the strategy for developing the 5CB CG force field. Notice that, to construct the optimal CG force field for 5CB liquid crystal molecule, we combined the thermodynamic quantities-based and structure-based CG methods and devised a careful guideline to parametrize the CG interaction potentials in a systematic way, and at same time attempted to use several fragment molecular systems to derive the CG nonbonded interaction parameters. Precisely for these reasons, the derived CG potentials could not just properly maintain the 5CB CG model thermodynamically and structurally consistent with the underlying atomistic model, but more importantly be transferable to the chemical-structure-like 6CB and 8CB liquid crystal systems across different thermodynamic conditions we studied. In particular, like the underlying atomistic models, the most favorable geometrical arrangement between the nearest neighboring 6CB or 8CB molecules in the nematic or/and smectic phases is also antiparallel in the CG simulations. Thus, although the dipolar interactions are not considered explicitly in the CG models, our strategy for constructing the CG force field makes the local structure information associated with the electrostatic interactions incorporating into the 5CB CG force field to some extent. Meanwhile, our results also show the limitation of the transferability of current 5CB CG force field. We found that, despite the formation of a smectic 8CB phase being seen in the CG simulations, the phase transition from nematic to smectic is postponed to the lower temperature at which CG 8CB systems

Table 3. Diffusion Coefficients of 8CB Systems from CG Simulations, Atomistic Simulations, and Experimental Measurements translatioal diffusion coefficients

T (K)

CG simulation

atomistic simulation

Diso (10−11 m2/s) D∥ (10−11 m2/s) D⊥ (10−11 m2/s) D∥/D⊥ D∥ (10−11 m2/s) D⊥ (10−11 m2/s) D∥/D⊥

330 305 305 305 280 280 280

18.2 4.5 1.5 3.0 0.7 0.4 1.75

16.3 15.1 7.8 1.93 1.9 1.1 1.73

experiment64 8.0 5.5 (310 K) 3.0 (310 K) 1.83 (310 K)

to the case of 6CB, the dynamics speedup also does not show up in the CG simulation of 8CB. We see that, for an isotropic phase at 330 K, the resulting Diso values from CG and atomistic simulations are 18.2 × 10−11 and 16.3 × 10−11 m2/s, respectively, which fall well within the range of corresponding experimental data. For a nematic phase at 305 K, we find that the CG simulated value D∥ (parallel to the director) is 4.5 × 10−11 m2/s, which is close to the experimental value of 5.5 × 10−11 m2/s at a nearby temperature of 310 K and the corresponding atomistic simulation result of 15.1 × 10−11 m2/s, while the CG simulated D⊥ (perpendicular to the director) is 1.5 × 10−11 m2/s, which is also in the correct range of the experimental measure at 310 K and atomistic simulated one at 310 K. As to the ratio of D∥/D⊥, the corresponding value is 3.0 from CG simulation, which can be well compared with the experimental and atomistic simulated values of 1.83 and 1.93. At transition into the smectic phase at 280 K, we observe a drop in both CG simulated D∥ and D⊥ and their values are well within the range of corresponding atomistic simulated results for a smectic phase at 280 K. As discussed above for 6CB, the 4657



Article

thermodynamic properties) parametrization route together with fragment-based method maximizes the possibility to simultaneously reproduce multiple physical properties (e.g., the phase behavior, some structural properties as well as diffusion coefficient ratios) of a given molecule or related molecules with similar chemical structures at other state points.

have higher densities and higher effective length-to-width ratios, and the resulting smectic phase structure is single-layer-like rather than partially interdigitated bilayer-like as observed in underlying atomistic model. These results remind us that we should pay more attention when applying the CG force field in the state points (e.g., smectic phase) significantly different from which the force field is explicitly parametrized for (e.g., nematic phase). In principle, the procedure to split the target molecule into small fragments (building blocks) and derive the CG nonbonded interaction parameters from isotropic liquids of these fragmental molecules is expected to maintain the transferability of CG parameters to some extent. According to our coarse-graining technique, among various CG bead pairs, only the nonbonded parameters of P3−P3 pair and C6−C6 pair are optimized in order to properly match a set of quantities such as bulk density, radial distribution functions, and nematic phase of 5CB systems at 300 K and 1 atm. Nevertheless, on the basis of such a scheme in the CG parametrization process, the current 5CB CG force field could inevitably incorporate the detailed structural information regarding the specific nematic phase of 5CB at 300 K and 1 atm, which may not allow for a locally very tight packing of 8CB molecules and thus fails to reproduce the smectic phase at 300 K and 1 atm as the underlying atomistic model. In general, increasing density can enhance the local packing of the molecules, and thus, it is not surprising to see TSN for the CG model at lower temperatures than for the chemically realistic one. Additionally, some necessary simplifications or approximations inherent in the CG process, such as the use of soft nonbonded potential and the removal of intramolecular barriers against rotation, which are expected to lower the effective repulsion between CG beads and reduce the shape anisotropy of the biphenyl segment, are also responsible for this fail of formation of the smectic phase at 300 K and 1 atm. Given this, increasing density or lowering temperature can strengthen excluded volume interactions and enhance the anisotropy of molecular shape, and the formation of the smectic phase is therefore pushed to the high density or low temperature end of the phase diagram when using current 5CB CG force field. Apart from these, electrostatic interaction, which is not explicitly modeled, may also limit the scope of our 5CB CG force field to precisely reproduce the experimental TSN as the underlying atomistic model. Similarly, the inherent simplifications and necessary scheme in the creation process of CG force field for example, choosing symmetric CG potentials which do not explicitly include electrostatic interactions and are parametrized by reproducing the target bulk density and RDF of the nematic 5CB phase at 300 K and 1 atm, as well as using soft nonbonded potential and excluding torsion barriers which will more or less increase the molecular flexibility of CG model and facilitate the interpenetration of smectic layers, may make it impossible for the CG 8CB model to precisely reproduce the bilayer-like smectic structures as the underlying atomistic model. Interestingly, in spite of this CG force field incorporating both thermodynamic and local structural information on the nematic 5CB phase at 300 K, it is transferable to study the dynamic properties of nCB isotropic, nematic, and/or smectic phases. For example, similar to the case of CG 5CB simulations, the anisotropic diffusion coefficient ratios for the isotropic, nematic, or/and smectic phases of both 6CB and 8CB systems were in the correct range of the experimental measurements and atomistic simulation results. These findings further suggest that the multiproperty (e.g., several key structural and

■

ASSOCIATED CONTENT

* Supporting Information S

RDFs for the pairs of CG sites of C3−C3, P3−P3, and C6−C6 obtained from CG and atomistic MD simulations of 6CB systems at T = 300 and 305 K, and RDFs for the pairs of CG sites of CN−CN, C2−C2, C3−C3, P3−P3, and C6−C6 obtained from CG and atomistic MD simulations of 8CB systems at T = 280 K (smectic phase) and 305 K (nematic phase). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 86-10-82618124. Fax: 8610-62559373. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful for the support of the NSF of China (20874110 and 20674093). We also thank the supercomputing center of CAS for supporting computing resources.

■

REFERENCES

(1) Zhang, J.; Su, J.; Ma, Y.; Guo, H. Coarse-Grained Molecular Dynamics Simulations of the Phase Behavior of the 4-Cyano-4′pentylbiphenyl Liquid Crystal System. J. Phys. Chem. B 2012, 116, 2075−2089. (2) Izvekov, S.; Voth, G. A. A multiscale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109, 2469−2473. (3) Srinivas, G.; Klein, M. L. Molecular Dynamics Simulations of Self-Assembly and Nanotube Formation by Amphiphilic Molecules in Aqueous Solution: a Coarse-Grain Approach. Nanotechnology 2007, 18, 205703. (4) Shelley, J. C.; Shelley, M. Y.; Reeder, R. C.; Bandyopadhyay, S.; Moore, P. B.; Klein, M. L. Simulations of Phospholipids Using a Coarse Grain Model. J. Phys. Chem. B 2001, 105, 9785−9792. (5) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (6) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K. Hierarchical Modeling of Polystyrene: From Atomistic to Coarse-Grained Simulations. Macromolecules 2006, 39, 6708−6719. (7) Harmandaris, V. A.; Reith, D.; Van der Vegt, N. F. A.; Kremer, K. Comparison between Coarse-Graining Models for Polymer Systems: Two Mapping Schemes for Polystyrene. Macromol. Chem. Phys. 2007, 208, 2109−2120. (8) Reith, D.; Putz, M.; Muller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624−1636. (9) Srinivas, G.; Discher, D. E.; Klein, M. L. Self-Assembly and Properties of Diblock Copolymers by Coarse-Grain Molecular Mynamics. Nat. Mater. 2004, 3, 638−644. (10) Milano, G.; Goudeau, S.; Müller-Plathe, F. Multicentered Gaussian-based Potentials for Coarse-Grained Polymer Simulations: Linking Atomistic and Mesoscopic Scales. J. Polym. Sci., Part B: Polym. Phys. 2005, 43, 871−885.

4658



Article

(11) Milano, G.; Müller-Plathe, F. Mapping Atomistic Simulations to Mesoscopic Models: A Systematic Coarse-Graining Procedure for Vinyl Polymer Chains. J. Phys. Chem. B 2005, 109, 18609−18619. (12) Li, X.; Ma, X.; Huang, L.; Liang, H. Developing Coarse-Grained Force Fields for cis-Poly(1,4-butadiene) from the Atomistic Simulation. Polymer 2005, 46, 6507−6512. (13) Li, X.; Kou, D.; Rao, S.; Liang, H. Developing a Coarse-Grained Force Field for the Diblock Copolymer Poly(styrene-b-butadiene) from Atomistic Simulation. J. Phys. Chem. B 2006, 124, 204909. (14) Müller-Plathe, F. Coarse-Graining in Polymer Simulation: From the Atomistic to the Mesoscopic Scale and Back. ChemPhysChem 2002, 3, 754−769. (15) Lyubartsev, A. P.; Laaksonen, A. Calculation of Effective Interaction Potentials from Radial-Distribution Functions - a Reverse Monte-Carlo Approach. Phys. Rev. E 1995, 52, 3730−3737. (16) Peter, C.; Kremer, K. Multiscale Simulation of Soft Matter Systems - from the Atomistic to the Coarse-Grained Level and Back. Soft Matter 2009, 5, 4357−4366. (17) Megariotis, G.; Vyrkou, A.; Leygue, A.; Theodorou, D. N. Systematic Coarse Graining of 4-Cyano-4 ′-pentylbiphenyl. Ind. Eng. Chem. Res. 2011, 50, 546−556. (18) Mukherjee, B.; Delle Site, L.; Kremer, K.; Peter, C. Derivation of Coarse Grained Models for Multiscale Simulation of Liquid Crystalline Phase Transitions. J. Phys. Chem. B 2012, 116, 8474−8484. (19) Izvekov, S.; Voth, G. A. Multiscale Coarse Graining of LiquidState Systems. J. Chem. Phys. 2005, 123, 134105. (20) Noid, W. G.; Chu, J.-W.; Ayton, G. S.; Krishna, V.; Izvekov, S.; Voth, G. A.; Das, A.; Andersen, H. C. The Multiscale Coarse-Graining Method. I. A Rigorous Bridge between Atomistic and Coarse-Grained Models. J. Chem. Phys. 2008, 128, 244114. (21) Noid, W. G.; Liu, P.; Wang, Y.; Chu, J.-W.; Ayton, G. S.; Izvekov, S.; Andersen, H. C.; Voth, G. A. The Multiscale CoarseGraining Method. II. Numerical Implementation for Coarse-Grained Molecular Models. J. Chem. Phys. 2008, 128, 244115. (22) Krishna, V.; Noid, W. G.; Voth, G. A. The Multiscale CoarseGraining Method. IV. Transferring Coarse-Grained Potentials between Temperatures. J. Chem. Phys. 2009, 131, 024103. (23) Liu, P.; Izvekov, S.; Voth, G. A. Multiscale Coarse-Graining of Monosaccharides. J. Phys. Chem. B 2007, 111, 11566−11575. (24) Wang, Y.; Feng, S.; Voth, G. A. Transferable Coarse-Grained Models for Ionic Liquids. J. Chem. Theory Comput. 2009, 5, 1091− 1098. (25) Hills, R. D.; Lu, L.; Voth, G. A. Multiscale Coarse-Graining of the Protein Energy Landscape. PLoS Comput. Biol. 2010, 6, e1000827. (26) Srinivas, G.; Klein, M. L. Coarse-Grain Molecular Dynamics Simulations of Diblock Copolymer Surfactants Interacting with a Lipid Bilayer. Mol. Phys. 2004, 102, 883−889. (27) Lopez, C. A.; Rzepiela, A. J.; de Vries, A. H.; Dijkhuizen, L.; Hunenberger, P. H.; Marrink, S. J. Martini Coarse-Grained Force Field: Extension to Carbohydrates. J. Chem. Theory Comput. 2009, 5, 3195−3210. (28) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S.-J. The MARTINI Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819− 834. (29) Peter, C.; Kremer, K. Multiscale Simulation of Soft Matter Systems. Faraday Discuss. 2009, 144, 9−24. (30) Engin, O.; Villa, A.; Peter, C.; Sayar, M. A Challenge for Peptide Coarse Graining: Transferability of Fragment-Based Models. Macromol. Theor. Simul. 2011, 20, 451−465. (31) Louis, A. Beware of Density Dependent Pair Potentials. J. Phys.: Condens. Matter 2002, 14, 9187. (32) Carbone, P.; Varzaneh, H. A. K.; Chen, X.; Müller-Plathe, F. Transferability of coarse-grained force fields: The Polymer Case. J. Chem. Phys. 2008, 128, 064904. (33) Qian, H.-J.; Carbone, P.; Chen, X.; Karimi-Varzaneh, H. A.; Liew, C. C.; Müller-Plathe, F. Temperature-Transferable CoarseGrained Potentials for Ethylbenzene, Polystyrene, and Their mixtures. Macromolecules 2008, 41, 9919−9929.

(34) Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; AlaNissila, T. Coarse-Graining Polymers with the MARTINI Force-Field: Polystyrene as a Benchmark Case. Soft Matter 2011, 7, 698−708. (35) Allen, E. C.; Rutledge, G. C. A Novel Algorithm for Creating Coarse-Grained, Density Dependent Implicit Solvent Models. J. Chem. Phys. 2008, 128, 154115. (36) Allen, E. C.; Rutledge, G. C. Coarse-Grained, Density Dependent Implicit Solvent Model Reliably Reproduces Behavior of a Model Surfactant System. J. Chem. Phys. 2009, 130, 204903. (37) Nielsen, S. O.; Lopez, C. F.; Srinivas, G.; Klein, M. L. A Coarse Grain Model for n-Alkanes Parameterized from Surface Tension Data. J. Chem. Phys. 2003, 119, 7043−7049. (38) Lorentz, H. A. Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Ann. Phys. 1881, 12, 127−136. (39) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (40) Nose, S. A Molecular-Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (41) Parrinello, M.; Rahman, A. Polymorphic Transitions in SingleCrystals - a New Molecular-Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (42) Luckhurst, G. R.; Veracini, C. A. E., Eds. The Molecular Dynamics of Liquid Crystals; Kluwer Academic Publishers: Dordrect, The Netherlands, 1994; Vol. 431. (43) Beguin, A. Sources of Thermodynamic Data on Mesogens; Gordon and Breach Science Publishers: New York, 1984. (44) deGennes, P. G.; Prost, J. The Physics of Liquid Crystals, second ed.; Oxford University Press: Oxford, U.K., 1993. (45) Kobinata, S.; Kobayashi, T.; Yoshida, H.; Chandani, A. D. L.; Maeda, S. Molecular-Conformation and Orientational Order in nCB Liquid-Crystals. J. Mol. Struct. 1986, 146, 373−382. (46) Beguin, A.; Dubois, J. C.; Lebarny, P.; Bonamy, F.; Buisine, J. M.; Cuvelier, P. Sources of Thermodynamic Data on Mesogens. Mol. Cryst. Liq. Cryst. 1984, 115, 119. (47) Miguel, E. d.; Rio, E. M. d.; Brown, J. T.; Allen, M. P. Effect of the Attractive Interactions on the Phase Behavior of the Gay−Berne Liquid Crystal Model. J. Chem. Phys. 1996, 105, 4234−4249. (48) Cifelli, M.; De Gaetani, L.; Prampolini, G.; Tani, A. Atomistic Computer Simulation and Experimental Study on the Dynamics of the n-Cyanobiphenyls Mesogenic Series. J. Phys. Chem. B 2008, 112, 9777−9786. (49) Izvekov, S.; Voth, G. A. Modeling Real Dynamics in the CoarseGrained Representation of Condensed Phase Systems. J. Chem. Phys. 2006, 125, 151101. (50) Lyubimov, I.; Guenza, M. G. First-Principle Approach to Rescale the Dynamics of Simulated Coarse-Grained Macromolecular liquids. Phys. Rev. E 2011, 84, 031801. (51) Matysiak, S.; Clementi, C.; Praprotnik, M.; Kremer, K.; Delle Site, L. Modeling Diffusive Dynamics in Adaptive Resolution Simulation of Liquid Water. J. Chem. Phys. 2008, 128, 024503. (52) Fu, C.-C.; Kulkarni, P. M.; Shell, M. S.; Leal, L. G. A Test of Systematic Coarse-Graining of Molecular Dynamics Simulations: Transport Properties. J. Chem. Phys. 2013, 139, 094107. (53) Junghans, C.; Praprotnik, M.; Kremer, K. Transport Properties Controlled by a Thermostat: An Extended Dissipative Particle Dynamics Thermostat. Soft Matter 2008, 4, 156−161. (54) Leadbetter, A. J.; Durrant, J. L. A.; Rugman, M. Density of 4 Normal-Octyl-4′-Cyano-Biphenyl (8CB). Mol. Cryst. Liq. Cryst. 1977, 34, 231−235. (55) Hadjsahraoui, A.; Louis, G.; Mangeot, B.; Peretti, P.; Billard, J. Liquid-Crystal Phase-Transitions of Thin-Layers - a Photothermal Analysis. Phys. Rev. A 1991, 44, 5080−5086. (56) Zhang, J.; Su, J.; Guo, H. An Atomistic Simulation for 4-cyano4′-Pentylbiphenyl and its Homologue with a Reoptimized Force Field. J. Phys. Chem. B 2011, 115, 2214−2227. (57) Gray, G.-W.; Lydon, J. New Type of Smectic Mesophase? Nature 1974, 252, 221−222. 4659



Article

(58) Fritz, D.; Koschke, K.; Harmandaris, V. A.; van der Vegt, N. F.; Kremer, K. Multiscale Modeling of Soft Matter: Scaling of Dynamics. Phys. Chem. Chem. Phys. 2011, 13, 10412−10420. (59) Onsager, L. The Effects of Shape on the Interaction of Colloidal Particles. Ann. N.Y. Acad. Sci. 1949, 51, 627−659. (60) Zhang, Z.; Guo, H. The Phase Behavior, Structure, and Dynamics of Rodlike Mesogens with Various Flexibility Using Dissipative Particle Dynamics Simulation. J. Chem. Phys. 2010, 133, 144911. (61) Williamson, D. C.; Jackson, G. Liquid Crystalline Phase Behavior in Systems of Hard-Sphere Chains. J. Chem. Phys. 1998, 108, 10294. (62) De Gaetani, L.; Prampolini, G. Computational Study Through Atomistic Potentials of a Partial Bilayer Liquid Crystal: Structure and Dynamics. Soft Matter 2009, 5, 3517−3526. (63) Wilson, M. Molecular Dynamics Simulation of Semi-flexible Mesogens. Mol. Phys. 1994, 81, 675−690. (64) Dvinskikh, S. V.; Furo, I.; Zimmermann, H.; Maliniak, A. Anisotropic Self-Diffusion in Thermotropic Liquid Crystals Studied by H-1 and H-2 Pulse-Field-Gradient Spin-echo NMR. Phys. Rev. E 2002, 65, 061701.

4660


Transferability of Coarse-Grained Force Field for nCB Liquid Crystal

Recommend Documents