Martin
Walker
and
Mark R.
Wilson
*
Durham University, Department of Chemistry, Lower Mountjoy, South Road, Durham DH1 3LE, UK. E-mail: mark.wilson@durham.ac.uk; Tel: +44 (0) 191 384 1158
First published on 29th September 2016
A simple dissipative particle dynamics (DPD) model is introduced, which can be used to represent a broad range of calamitic mesogens. The model allows for antiparallel association that occurs naturally in a number of mesogens with terminal dipoles, including the 4-n-alkyl-4′-cyanobiphenyl (nCB) series. Favourable antiparallel interactions lead to the formation of SmAd phases in which the layer spacing is intermediate between monolayer and bilayer. The model is easily tuned to vary the strength of antiparallel association and the SmA layer spacing, and to give either isotropic–smectic or isotropic–nematic–smectic phase sequences. The model allows for a range of other smectics: including SmA1 phases exhibiting microphase separation within layers, and smectics A structures with more complicated repeat units. For large system sizes (≥50000 molecules) in the nematic phase, we are able to demonstrate the formation of three distinct types of cybotactic domains depending on the local interactions. Cybotactic domains are found to grow in the nematic–smectic pretransitional region as the system moves closer to TSN.
Molecular simulations have proven to be critical in explaining some of the fundamental aspects of liquid crystal phase behaviour.8,9 Early simulation studies used simple hard repulsive shapes to model liquid crystal molecules, e.g. spherocylinders to represent prolate molecules10,11 and cut disks to represent oblate molecules.12 Hard potentials were able to generate simple phase diagrams with isotropic, nematic and smectic phases, employing only an entropic contribution to the total free energy. As computational power increased the development of anisotropic-attractive-coarse-grained models, such as the Gay-Berne potential,13–15 led to increased use of simulations to explore the link between molecular structure, phase behaviour and the bulk properties of liquid crystal phases.16–18 For nematics, it is now possible to generate very good predictions for transition temperatures from atomistic simulations of liquid crystal phases.19–21 However, for more complex mesophases, coarse-grained models are extremely useful; particularly because the system sizes and time scales required are usually too demanding for atomistic level models.
While excluded volume factors play a major role in driving mesophase formation, specific pair-wise interactions have long been known to provide an important contribution to the formation of some liquid crystal phases.22–24 Moreover, subtle effects such as changes in the tilt angle within smectic C liquid crystals, and re-entrant phase behaviour can be attributed to specific (often dipole–dipole coupling) interactions.25
The layer structure of a number of smectic liquid crystals is also influenced by dipole–dipole coupling interactions. For example, 4-n-octyl-4′-cyano-biphenyl, 8CB, which has a large dipole at the end of the molecule, adopts a smectic-A phase (classified as SmAd), in which the layer spacing, dlayer = 1.4 molecular lengths.26–28 The extended layer spacing is attributed to a preferred antiparallel packing of molecules,29,30 predominately driven by antiparallel dipole–dipole ordering. Atomistic simulations have been able to reproduce both the dipole–dipole ordering and the layer spacing in 8CB,20,31,32 albeit for a small number of smectic layers (on account of system size limitations). The simulations show that the layer density wave is considerably broader than that observed in other smectic A phases. Moreover, despite the fact that other parts of the molecular interaction have a far stronger influence on transition temperatures than dipole–dipole interactions, the latter are vital in determining the overall structure of the phase.
The SmAd phase of 8CB is only one of the possible smectic A phases that can occur when there is competition between two competing incommensurate lengths: that of a molecule and that of a pair of molecules. de Gennes and Prost show, with a Landau theory of frustrated smectics,33 how this can lead to a series of different smectics A phases with different density waves: SmA1 (monolayer), SmA2 (bilayer), SmAd and SmÃ.
In this paper, we present a coarse-grained dissipative particle dynamics (DPD) model for liquid crystal molecules, in which the effects of specific molecular interactions (such as those leading to dipole–dipole correlation or to microphase separation in smectics) can be incorporated easily. We show that we can represent nematic and smectic SmAd phases within the model (mimicking the behaviour of 8CB); and we show how a wide range of smectic-A phases can be generated by tuning the local interaction potential to favour specific interactions leading to local (dynamic) molecular recognition. We show also that DPD (because of the large system sizes and time scales possible) provides a valuable tool for looking at pretransitional fluctuations and cybotactic behaviour within nematic systems.
In DPD three simple forces act on the system: a conservative force, FC, a dissipative force, FD, and a random force, FR.
(1) |
(2) |
(3) |
To provide a generic model for calamitic liquid crystals, (i.e. instead of a fully chemically tractable coarse graining scheme) we represent a liquid crystal mesogen by eight beads that are linearly bonded with harmonic bonds to maintain connectivity. Here,
(4) |
(5) |
Equations of motion were solved using the Velocity Verlet integration algorithm, using the DPD module of the DL_MESO50,51 program (DL_MESO_DPD) in the constant-NVT (canonical) ensemble. In each case, simulations were started from a random configuration (position and orientation) of 5000 molecules, which were initially equilibrated to an isotropic phase (T* = 2.0) and then cooled in a stepwise manner in increments of 0.1ΔT*/500000δt. We found that 500000 steps with a time step, δt, of at each temperature was sufficient to ensure a well equilibrated phase in all cases studied. For some simulations we used larger system sizes of 50000 molecules to check the structure of smectic phases with a larger number of layers. The DPD parameter for noise amplitude was taken to be σ = 3.6710. The temperature is controlled by γ, such that the relationship between γ and σ is maintained, (σ2 = 2γT*). All simulations were performed at a particle site density ρ = 3.0(l−3). A mesogen showing a clear isotropic–nematic–smectic transition sequence was found for a model of eight linearly bonded beads, of a single type (type A), and with conservative force parameter aAA = 45.0ε. We define this as model 0 and this acts as a basis of comparison for later smectic phases.
Normally, for simulations of layered liquid crystals with relative small systems sizes, constant-NpT simulations are often needed in order to avoid the formation of layer structures that are incommensurate with the box dimensions. However, for the DPD method we employ here this is not necessary. Firstly, we use large system sizes, where we can check that the simulation box does not impose structure on the phase. Secondly, the simulation times used are extremely long relative to standard molecular dynamics models (in terms of both molecular diffusion and reorientation), so the director is quite easily able to rotate to a different orientation should a new phase start to form that is initially incommensurate with the periodic boundary conditions.
To favour antiparallel molecular interactions, two additional bead types were introduced (types B (blue) and C (red)). For all models the interaction parameters for the new beads were set to aAA = aAB = aAC = aBB = aCC = 45.0ε. However, the cross term aBC used to favour antiparallel molecular recognition, was set at a reduced value (with values of between 35.0ε and 5.0ε used below). Hence the magnitude of aBC controls the preference for antiparallel interactions, and the position of B and C groups within the molecule changes the location of the interaction (mimicking dipolar regions in different parts of the molecule).
Thirteen separate models were studied, as shown schematically in Fig. 1. As the strength of preferred B–C interactions is always relative to the maximum value of aij(45.0ε), we can define a “molecular recognition” variable arec = aij(max) − aBC to measure this, and classify the region between and including B and C in a molecule as the “molecular recognition” site.
We note that in this simple model, molecular recognition is limited to interactions arising from just two sites per molecule, rather than from those arising from a large molecular segment. However, this is sufficient to enable favourable antiparallel configurations (in 8CB arising from π–π interactions in addition to antiparallel dipole interactions) without the use of computationally costly coulombic interactions.
The nematic order parameter for these systems was obtained by first defining a long-axis unit vector for each molecule, ui (from diagonalization of the molecular inertia tensor) followed by diagonalization of the order tensor
(6) |
rm‖,c‖ = rij·. | (7) |
Fig. 3 shows the parallel pair distribution functions g(rc‖), from which a single layer spacing of d = 4.15l, is clearly observable. At low temperatures (T* = 0.8) the mean bond length is equal to 0.52l. Hence d corresponds to approximately one molecular length and the phase is a standard monolayer smectic A phase (SmA1).
Using the mean bond length (above) together with a mean width of 0.8l (minimum width = 0.5l), we obtain a molecular aspect ratio of ≈5.6:1 for our reference model (slightly less than this value if the flexibility of bond angles is taken into account). The low temperature model system is therefore slightly shorter than a 5:1 spherocylinder model.11
The designation for the smectic phase is SmAd, and the model therefore mimics the smectic phase behaviour observed in longer chain homologues of the nCB (4-n-alkyl-4′-cyanobiphenyl) series, which show antiparallel dipole correlation arising from the large dipole associated with the terminal cyano group. This behaviour has been observed in atomistic simulation studies of 8CB (which has a layer spacing of 1.4 molecular lengths).20,31,32 However, the atomistic simulations are at the limit of what is currently computationally feasible. So for example, in ref. 32, pair distribution plots were only extendible to ∼30 Å, when the experimental layer spacing is 30.5 Å. Here, for simulations of 50000 molecules, we can resolve 4 layers in g(rm‖).
As in 8CB, the SmAd phase observed is not simply built from antiparallel dimers, in which every B site has two interacting C neighbours (in 2 dimensions). Such a phase would result in a shorter interlayer spacing (comparable to that of model 0) but with very tight (entropically unfavourable) packing of molecules from adjacent layers into spaces between dimers. Instead an alternative molecular packing is possible, as shown in the inset to Fig. 3, where half the C sites are able to interact with multiple B sites. This packing behaviour explains the additional 0.9l peak and the inter-layer spacing of 6.3l. It also leads to broader smectic layers than those seen in model 0.
Decreasing the molecular recognition parameter for model 1 to arec = 30.0ε, weakens the antiparallel association, resulting in subtly different phase behaviour. On cooling from the isotropic phase, a stable nematic phase is formed at T* = 1.1, prior to the formation of a SmAd phase. The local packing differs from that of a standard nematic, in that short range positional order is clearly detectable in g(rm‖), as shown in Fig. 4. We attribute the local positional order to transient antiparallel association, in agreement with the pre-transitional behaviour as noted in dielectric studies of cyanobiphenyl nematics.60–63 As temperature is reduced (Fig. 4), the strength of antiparallel correlation increases, resulting in cybotactic domains of SmAd phase within the nematic. Further reduction in the molecular recognition parameter arec = 10.0ε results in a nematic phase (Fig. 2), with far weaker local antiparallel correlation. At this value of arec the long layer spacing of a SmAd phase is lost, replaced by a standard SmA1 phase (at T* < 0.9), as observed in model 0.
Fig. 4 g(rm‖) for model 1 with arec = 30ε showing the growth of pretransitional fluctuations on cooling the nematic phase: T* = 1.1 (bold line), T* = 1.0 (dashed line) and T* = 0.9 (dotted line). |
Model 4 (Fig. 1), has an increased distance between the antiparallel sites, leading to subtly different phase behaviour in comparison to model 1. For 10ε ≤ arec ≤ 40ε, on cooling from the isotropic phase, a nematic with local antiparallel correlation is formed with an onset temperature of T* = 1.1, and a smectic phase is formed at T* = 1.0. Both g(rm‖) and g(rc‖) show peaks at both half integral and integral molecular lengths. The shorter correlation length can be attributed to strong antiparallel packing causing intercalation of adjacent layers. The inset in Fig. 5 shows the local packing within the model, which supports each B interacting with C sites, and vice versa. However, the strict repeat distance for layers along the director is equivalent to two molecular lengths. The half molecular length displacement can occur in either a +ve or −ve direction along the director, and so packing of molecules into layers gives relatively sharp intensity peaks in g(rm‖) in comparison to the diffuse peaks observed for model 1. It is important to note that unlike 8CB and model 1, where only two sub-layers contribute to a single full layer; model 4 has three contributing sub-layers, as demonstrated by a comparison of the schematic insets for Fig. 3 and 5.
The phase behaviour exhibited by model 4 is replicated by model 7 with the sole exception that for arec = 10.0ε the antiparallel correlation is sufficiently weak that it is lost completely in the smectic phase, with the formation of a SmA1 phase (as seen in model 0).
Model 9 exhibits two types of SmAd phase. For arec = 40.0ε and arec = 30.0ε, a strongly intercalated smectic phase is formed, comparable to that of model 4. However, for arec < 20.0ε a SmA1 phase is formed. The intercalated smectic phase has similarities with both models 1 and 4. Two peaks are seen in g(rm‖) at rm‖ = 2.7l and rm‖ = 5.2l, which correspond to the half layer and full layer peaks of model 4. These peaks are at slightly larger values of rm‖ than those observed in model 4, giving smectic layers with a width of 1.3 molecular lengths. The antiparallel packing is shown schematically in Fig. 6 and highlights that the packing arrangement adopted shares similarities with that of model 4, i.e. the B and C sites interact in an identical arrangement. However, the additional spacing between the two interaction sites subtly shifts the layer spacing to greater distances in a similar manner to model 1. Also the satellite peak, which occurred in model 1 at r = 0.9l, is displaced in model 9 towards the middle of the layer.
A summary of the phase behaviour of each of these models is given in Table 1.
Model number | a rec/ε | Mesophases | T NI | T SN | Smectic phase type | Main peaksa in g(rm‖)/molecular length |
---|---|---|---|---|---|---|
a Here, the main peaks in g(rm‖)/are measured for distances as defined in Fig. 1. They therefore do not necessarily correspond to observable X-ray peaks, given the anti-parallel orientation of some neighbouring molecules. | ||||||
0 | 0 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1 |
1 | 40 | iso–sm | 1.15 | SmAd | 1.5 | |
1 | 30 | iso–nem–sm | 1.15 | 0.95 | SmAd | 1.5 |
1 | 20–10 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1 |
2 | 40–10 | iso–nem–sm–sm | 1.05 | 0.95 | SmA1 | 1 |
3 | 40–10 | iso–nem–sm | 1.15 | 1.05 | SmA1mic | 1 |
4 | 40–10 | iso–nem–sm | 1.15 | 1.05 | SmA2 | 1, 0.5 |
5 | 40–10 | iso–sm–sm | 1.05 | SmA1 | 1 | |
5 | 30–10 | iso–nem–sm–sm | 1.15 | 0.95 | SmA1 | 1 |
6 | 40–30 | iso–sm | 1.05 | SmA1mic | 1 | |
6 | 20–10 | iso–nem–sm | 1.05 | 0.95 | SmA1mic | 1 |
7 | 40–20 | iso–nem–sm | 1.15 | 0.95 | SmA2 | 1, 0.5 |
7 | 10 | iso–nem–sm | 1.15 | 0.85 | SmA1 | 1, 0.5 |
8 | 40 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1, 0.5 |
8 | 30–10 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1 |
9 | 40–30 | iso–sm | 1.15 | SmAd | 1.3, 0.7 | |
9 | 30–10 | iso–nem–sm | 1.15 | 0.95 | SmAd | 1.3, 0.7 |
9 | 20–10 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1 |
10 | 40–30 | iso–sm | 1.05 | SmA1mic | 1 | |
10 | 20–10 | iso–sm–nem | 1.15 | 0.95 | SmA1mic | 1 |
11 | 40–10 | iso–nem–sm | 1.15 | 0.95 | SmA1 | 1 |
12 | 40–10 | iso–sm | 1.55 | SmA1mic | 1 |
Fig. 7 Left: g(rm‖) for model 6 with arec = 20ε, smectic phase at T* = 0.9l (bold line), nematic phase at T* = 1.0l (dotted line). Right: A simulation snapshot of the smectic phase for model 6. |
Model 10 behaves in an almost identical way to model 6, but with the position of the microphase separated domains shifted slightly. Model 12 behaves in a similar way, but with the positions of the B and C sites at the end of the molecule providing an additional inter-layer interaction. This leads to an exceptionally stable smectic A1 phase (stable up to T* ∼ 1.55), with no nematic phase observable at any temperatures for 10ε ≤ arec ≤ 40ε.
Models 3 and 11 differ slightly from the others in this section because the B and C sites are not equidistant from the centre of the molecule. This loss of symmetry, causes the smectic layers to be less strongly ordered (relative to models 6 and 12 respectively), reducing the thermal stability of the SmA1 phase relative to the nematic (see Table 1).
The two smectic phases are easily observed in model 2 (with arec = 40ε), which exhibits a smectic to smectic phase transition at T* ≈ 0.85. The change in structure of the smectic is shown in g(rm‖) plots (Fig. 8), where at the phase transition additional peaks appear. In the schematic diagram (right of Fig. 8) we show how the additional correlations, corresponding to these peaks, arise from the local packing arrangement of molecules. In the case of the lower temperature phase, there is no specific preference for an antiparallel ordering of neighbours; instead, additional correlations arise from both parallel and antiparallel neighbours in equal measure.
In the high-temperature intercalated smectic phase, the B and C sites do give rise to antiparallel ordering (see the top schematic diagram in Fig. 8), generating diffuse smectic layers. Here, however the average spacing remains equal to one molecular length. This complex layer structure is entropically favoured; hence the transition to a more conventional SmA1 phase as temperature is reduced.
The phase behaviour of model 5 is very similar to model 2 with small differences in transition temperatures and nematic stability. However, the behaviour of model 8 is subtly different. For arec = 40.0ε, isotropic–nematic (T* = 1.0)–smectic (T* = 0.9) transitions are seen, with the smectic phase exhibiting intercalated layers as observed in model 4. Here, the g(rm‖) function shows peaks corresponding to a half molecular length and a full molecular length (similar to model 4). The nematic phase that forms is very sensitive to arec. For arec ≥ 30ε, clear cybotactic domains are formed, which are ferroelectric in nature (see left hand side of the inset in Fig. 9). At lower values of arec cybotactic domains remain but the domains are of SmA1 character and are not ferroelectric (see right hand side of the inset in Fig. 9). The differences in the cybotactic domains show up as differences in the g(rm‖) plots in Fig. 9, noting that the peak intensities are much smaller than seen in smectic phases.
The model also demonstrates the formation of cybotactic nematics in which pretransitional smectic fluctuations occur. For different models we identified three distinct types of cybotactic domains:
• normal SmA1 regions,
• SmA regions with molecules showing strong antiparallel correlation,
• ferroelectric smectic regions.
As expected, cybotactic domains grow as the system moves closer to the smectic–nematic phase transition.
Recently, there has been considerable controversy arising from the role (or otherwise) of cybotactic domains within the nematic phase. In particular, for some bent-core nematics these domains appear to be very extensive and different domain structures are possible;64 and it has been suggested that cybotactic domains may account for the NMR biaxiality observed in nematic phases of some mesogens.65 This current work suggests that simple DPD models can represent the twin features of molecular shape and specific pair-wise interactions responsible for cybotactic domain formation. Work applying these types of models to bent-core mesogens is currently under way in our laboratory.
This journal is © The Royal Society of Chemistry 2016 |