Simulation insights into the role of antiparallel molecular association in the formation of smectic A phases

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 0 -cyanobiphenyl ( n CB) series. Favourable antiparallel interactions lead to the formation of SmA d 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 SmA 1 phases exhibiting microphase separation within layers, and smectics A structures with more complicated repeat units. For large system sizes ( Z 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 T SN .


Introduction
Thermotropic liquid crystals represent a fascinating area of soft matter science.High interest arises from the wide range of different liquid crystal mesophases that have been produced by Chemists; and also from the large number of commercial applications.The former includes new forms of nematic phase, 1 and a preponderance of smectic phases of varying degrees of translation and rotational order. 2 Applications cover the areas of thermochromic materials, 3 lasers, 4 high-tech lubricants, 5 liquid crystal displays, 6 and a range of optoelectronic applications such as optical filters and switches, beam-steering devices and spatial light modulators. 7olecular simulations have proven to be critical in explaining some of the fundamental aspects of liquid crystal phase behaviour. 8,9Early simulation studies used simple hard repulsive shapes to model liquid crystal molecules, e.g.spherocylinders to represent prolate molecules 10,11 and cut disks to represent oblate molecules. 12Hard potentials were able to generate simple phase diagrams with isotropic, nematic and smectic phases, employing only an entropic contribution to the total free energy.0][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.
3][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. 25he layer structure of a number of smectic liquid crystals is also influenced by dipole-dipole coupling interactions.7][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 SmA d 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: SmA 1 (monolayer), SmA 2 (bilayer), SmA d and SmA ˜.
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 SmA d 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.

Computational methods
Originally designed for block co-polymer systems, DPD has a smooth repulsive potential that allows for a relatively large time step.Coupling the reduced degrees of freedom of a coarse grained simulation with a large time step has enabled simulation of a wide range of soft matter: block co-polymers, [34][35][36] colloidal suspensions, 37 lipid bilayers, 38,39 nanoparticles, 40,41 and liquid crystalline phases, 42,43 where the hydrodynamic time-scale can become important.Molecular systems are readily coarse-grained to a DPD representation, replacing multiple atoms or chemical groups with single-site beads.5][46][47][48] Here, we simulate preferred antiparallel molecular interactions by using two sites, B and C, which have a preferred B-C interaction.
In DPD three simple forces act on the system: a conservative force, F C , a dissipative force, F D , and a random force, F R .
where a ij represents the maximum repulsive force between particles i and j, and r cut is the interaction cut-off, which, as standard, is taken to be equal to a single unit of length.
where o = r cut À |r ij | and the relationship between s and g (s 2 = 2gk B T = 2gT*) effectively controls the temperature.Equilibration rates with respect to temperature can be controlled with the friction coefficient o, while y is a random normally-distributed variable (0 o y o 1), enabling F R to mimic the influence of atomic motions within the larger bead representation.Further details of the technique have been reviewed in the literature.
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, where k bond is the force constant, expressed in DPD units (scaled mass, m = 1, distance expressed relative to the size of a DPD bead l = r cut = 1, time in reduced seconds s, and energies scaled by unit energy, e = 1) takes a value of 20el À2 with an equilibrium distance of zero.The repulsive component of the bond is imparted from the conservative force, F C .The shape of a mesogen is maintained with harmonic angle springs: where k angle , the angle force constant, takes the value k angle = 20(erad À2 ) and an equilibrium angle of p radians.Equations of motion were solved using the Velocity Verlet integration algorithm, using the DPD module of the DL_MESO 50,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 molecules, which were initially equilibrated to an isotropic phase (T* = 2.0) and then cooled in a stepwise manner in increments of 0.1DT*/500 000dt.We found that 500 000 steps with a time step, dt, of t 0 l=ðm=eÞ ¼ 0:02 at each temperature was sufficient to ensure a well equilibrated phase in all cases studied.For some simulations we used larger system sizes of 50 000 molecules to check the structure of smectic phases with a larger number of layers.The DPD parameter for noise amplitude was taken to be s = 3.6710.The temperature is controlled by g, such that the relationship between g and s is maintained, (s 2 = 2gT*).All simulations were performed at a particle site density r = 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 a AA = 45.0e.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 a AA = a AB = a AC = a BB = a CC = 45.0e.However, the cross term a BC used to favour antiparallel molecular recognition, was set at a reduced value (with values of between 35.0e and 5.0e used below).Hence the magnitude of a BC 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 a ij (45.0e), we can define a ''molecular recognition'' variable a rec = a ij (max) À a BC 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 p-p 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, u i (from diagonalization of the molecular inertia tensor) followed by diagonalization of the order tensor where (a, b = x, y, z), N is number of molecules and d ab is the Kronecker delta.The nematic order parameter, S 2 , was taken as the largest eigenvalue, with the director for the phase, n ˆ, given by the associated eigenvector.To further characterise mesophases, we introduce two pair correlation functions, g(r c 8 ) and g(r m 8 ), measured for intermolecular separations r c 8 and r m 8 along n ˆ.Here, where r ij is the intermolecular vector between molecules i and j measured for r c 8 relative to the geometric centres of molecules, and for r m 8 measured relative to the centres of the ''molecular recognition'' site, i.e. between beads B and C (as shown in Fig. 1).Noting that the latter is different for each individual model.Perpendicular radial distribution functions, g(r c > ), g(r m > ), were also calculated for distances perpendicular to the director.These allowed us to check the fluidity of smectic layers.For the systems studied here, we concentrate discussion on temperatures of T* Z 0.8, where g(r c > ) curves show systems to be fluid for each of the thirteen models and four values of a rec .Hence, we restrict ourselves to studying nematic and smectic-A phases, rather than more higher ordered smectics or crystals.

Results and discussion
3.1 Reference system -model 0 (SmA 1 phase) The reference model for this study, model 0, is provided by a mesogen of eight linearly-attached type A beads (Fig. 1).As expected, from previous work on hard sphere chains, [52][53][54] hard and soft spherocylinders, 11,[55][56][57][58][59] and DPD chains, 42 model 0 exhibits isotropic, nematic and smectic-A phases.On cooling from T* = 2.0 the isotropic-nematic transition was characterised by a discontinuous change in order parameter (see Fig. 2) at T* E 1.15, and the onset of the smectic phase (at T* E 0.9) was monitored through the growth of peaks in the radial distribution function resolved parallel to the director.Fig. 3 shows the parallel pair distribution functions g(r c 8 ), 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 (SmA 1 ).
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 E5.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. 112 Interdigitated bilayers (SmA d phases): models 1, 4, 7 and 9 The models of this class all have a ''molecular recognition'' region located at one end of the molecule.Model 1 shows the presence of different smectic phases, and also shows a high sensitivity to the magnitude of a rec .Fig. 2 shows the phase sequence obtained for four values of a rec .For the highest value considered (a rec = 40.0e),a single transition is observed, from isotropic to smectic at T* = 1.1.The smectic A phase formed differs from that of model 0, as the layer spacing is considerably larger than one molecular length.This is evident from the pair correlation plots shown in Fig. 3. Here, the choice of pair correlation function is clearly important.g(r c 8 ) provides minimal information but g(r m 8 ) shows a repeat layer spacing of E1.5 molecular lengths.The centres of the layers are provided by the ''molecular recognition'' region, with the initial peak in g(r m 8 ) found at r m 8 = 0, together with a shoulder at 0.9l.From the snapshot, also shown in Fig. 3, pair-wise antiparallel correlation can be clearly observed within microphase segregated regions.
The designation for the smectic phase is SmA d , and the model therefore mimics the smectic phase behaviour observed in longer chain homologues of the nCB (4-n-alkyl-4 0 -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,32However, 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 B30 Å, when the experimental layer spacing is 30.5 Å.Here, for simulations of 50 000 molecules, we can resolve 4 layers in g(r m 8 ).
As in 8CB, the SmA d 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 a rec = 30.0e,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 SmA d phase.1][62][63] As temperature is reduced (Fig. 4), the strength of antiparallel correlation increases, resulting in cybotactic domains of SmA d phase within the nematic.Further reduction in the molecular recognition parameter a rec = 10.0eresults in a nematic phase (Fig. 2), with far weaker local antiparallel correlation.At this value of a rec the long layer spacing of a SmA d phase is lost, replaced by a standard SmA 1 phase (at T* o 0.9), as observed in model 0.
Model 4 (Fig. 1), has an increased distance between the antiparallel sites, leading to subtly different phase behaviour in Fig. 2 The temperature dependence of the nematic order parameter, S 2 , for model 0 (black symbols), and model 1 with a rec = 40e (red symbols), a rec = 30e (green symbols), a rec = 20e (blue symbols), a rec = 10e (pink symbols).Squares show a stable SmA 1 phase, diamonds show a stable SmA d phase, triangles show a stable nematic phase and circles show a stable isotropic phase.comparison to model 1.For 10e r a rec r 40e, 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(r m 8 ) and g(r c 8 ) 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(r m 8 ) 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 a rec = 10.0e the antiparallel correlation is sufficiently weak that it is lost completely in the smectic phase, with the formation of a SmA 1 phase (as seen in model 0).
Model 9 exhibits two types of SmA d phase.For a rec = 40.0eand a rec = 30.0e,a strongly intercalated smectic phase is formed, comparable to that of model 4.However, for a rec o 20.0e a SmA 1 phase is formed.The intercalated smectic phase has similarities with both models 1 and 4. Two peaks are seen in g(r m 8 ) at r m 8 = 2.7l and r m 8 = 5.2l, which correspond to the half layer and full layer peaks of model 4.These peaks are at slightly larger values of r m 8 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.
3.3 SmA 1 phases with microphase separation and in-layer antiparallel ordering: models 3, 6, 10, 11 and 12 Models 3, 6, 10, 11 and 12 all follow very similar phase behaviour.In all these models SmA 1 phases are formed with a layer spacing of approximately one molecular length.Each system shows strong antiparallel ordering of molecules within the layer, and microphase separation into separate domains.Model 6 exhibits representative behaviour for this class of molecules.For a rec = 40.0eand a rec = 30.0ea single isotropic to SmA 1 (T* = 1.0) transition is observed upon cooling, with sharp peaks in g(r m 8 ) indicating that molecules are ''strongly bound'' to smectic layers and more restricted in their motion parallel to the director.This arises because of strong in-layer antiparallel correlation, which leads to five microphase separated regions A-BC-A-BC-A.This ordering occurs naturally when the molecules sit exactly side-by-side and head-to-tail (see the snapshot in Fig. 7).As the strength of the interaction between the two recognition sites is reduced (a rec = 20.0e) a nematic phase (T* B 1.05) is observable, with some (weak) positional order shown in g(r m 8 ), prior to smectic phase formation (T* B 0.95).
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 A 1 phase (stable up to T* B 1.55), with no nematic phase observable at any temperatures for 10e r a rec r 40e.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 SmA 1 phase relative to the nematic (see Table 1).

3.4
Complex SmA 1 phases without microphase separation: models 2, 5, 8 The final type of phase behaviour is covered by models 2, 5 and 8. Here, complex SmA 1 phases are formed with a layer spacing equivalent to one molecule.All these models form a standard SmA 1 without microphase separation.In addition, for higher values of a rec the models also exhibit a second smectic: in this case a microphase-separated intercalated smectic A phase, which retains a single layer repeat distance.
The two smectic phases are easily observed in model 2 (with a rec = 40e), which exhibits a smectic to smectic phase transition at T* E 0.85.The change in structure of the smectic is shown in g(r m 8 ) 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 a Here, the main peaks in g(r m 8 )/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.This journal is © The Royal Society of Chemistry 2016 molecular length.This complex layer structure is entropically favoured; hence the transition to a more conventional SmA 1 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 a rec = 40.0e,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(r m 8 ) 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 a rec .For a rec Z 30e, clear cybotactic domains are formed, which are ferroelectric in nature (see left hand side of the inset in Fig. 9).At lower values of a rec cybotactic domains remain but the domains are of SmA 1 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(r m 8 ) plots in Fig. 9, noting that the peak intensities are much smaller than seen in smectic phases.

Conclusions
We present a simple DPD model, which can be used to represent calamitic mesogens.The model allows for antiparallel association that occurs naturally in a wide range of mesogens with terminal dipoles, including the 4-n-alkyl-4 0 -cyanobiphenyl series.In particular the model allows for the simulation of SmA d phases in which the layer spacing is intermediate between monolayer and bilayer, with the stabilization of intermediate layer spacings controlled by the strength of antiparallel association.The phase sequences are easily controlled by a single parameter, a rec , which controls the strength of antiparallel association relative to kT.The model also allows for a range of other smectics: including SmA 1 phases exhibiting microphase separation within layers (and hence very stable smectics); and smectic A structures in which we see an overall repeat structure of two molecular lengths but with repeat distances in g(r m 8 ) of approximately half a molecular length.These are all obtainable by subtle changes in the position of the ''molecular recognition'' sites within the molecule.
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 SmA 1 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. 65This 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.

Fig. 1
Fig. 1 Left -Schematic diagram showing the thirteen DPD models studied in this work.Pink represents type A beads, blue type B and red type C. Beads are numbered sequentially from right to left.Right -Definition of parallel distances r c 8 and r m 8 for model 1.

Fig. 3
Fig. 3 Top left: g(r c 8 ) for the smectic-A phase of model 0 at T* = 1.1.Top right: A simulation snapshot showing the structure of the smectic-A phase of model 0 at T* = 1.1.Bottom left: g(r m 8 ) (black) and g(r c 8 ) (red) for the smectic-A d phase of model 1 at T* = 1.1, the solid vertical line indicates the expected spacing for a standard SmA 1 phase.A schematic showing the expected packing of molecules is shown in the inset.Bottom right: A simulation snapshot showing the structure of the smectic-A d phase of model 1 at T* = 1.1.

Fig. 5
Fig. 5 Left -g(r c 8 ) (red line) and g(r m 8 ) (black line) for model 4. Right -A simulation snapshot for model 4 with a rec = 30e at T* = 0.8.A schematic showing the expected packing of molecules is shown in the inset.

Fig. 6
Fig. 6 Left: g(r m 8 ) for model 9 with a rec = 30e (bold line) and a rec = 20e (dotted line).Right: Simulation snapshot for a rec = 30e.A schematic showing the expected packing of molecules for the a rec = 30e model is shown in the inset.

Fig. 8
Fig. 8 Left: g(r m 8 ) for model 2 with a rec = 40e, T* = 1.0l (dotted), T* = 0.9l (dashed), T* = 0.8l (bold).Right: Schematic showing the two packing arrangements, with indication of the expected g(r m 8 ) peaks.The vertical colour coded lines on the g(r m 8 ) functions indicate the expected location and intensity of the satellite (black) and full layer (blue) peaks.

Table 1
Summary showing the phase behaviour for all model studied.(Transition temperatures are given to T* values AE 0.05, ''mic'' indicates a smectic showing microphase separated domains)