Emergent chirality in achiral liquid crystals: insights from molecular simulation models of the behaviour of bent-core mesogens

ABSTRACT

We report insights into the behaviour of achiral bent-core liquid crystals using (a) atomistic models of the PnOPIMB series of mesogens and (b) simple coarse-grained models that demonstrate the spontaneous transfer of chirality between molecules by preferential selection of chiral conformations. The models explain the unusual phenomena that chiral interactions can lead to an increase in chirality (decrease in pitch) of a bulk cholesteric liquid crystal when doped with an achiral material. We demonstrate also that chiral interactions are sufficiently strong to induce nanophase separation of an achiral mesogen within a liquid crystal phase, into nanodomains of opposite handedness by preferential selection of conformations that are instantaneously chiral (i.e. nanophase separation of $R$ and $S$ conformations of a molecule that is strictly achiral). We suggest that this unusual behaviour contributes to the templating of chiral structures in bulk phases and, hence, to many of the extraordinary properties of bent-core mesogens that include the formations of helical nanofilaments and the chiral dark conglomerate phase.

Graphical Abstract

KEYWORDS:

• Chiral
• bent-core mesogens
• molecular modelling
• computer simulation

1. Prologue: a Festschrift for Claudio Zannoni

I was woken at 5.30am by gentle birdsong and the sun streaming through the wooden shutters of my small monastic cell. I awoke with the idea for this paper clearly in my mind: An idea born out of the conversations of the night before over Marsala and wonderful Italian hospitality.

Table 1. Optimised Ryckaert-Bellemans coefficients, $C_n$ $(\mathrm{k}\mathrm{J}\mathrm{m}\mathrm{o}{\mathrm{l}}^{-1})$, obtained from fitting to quantum chemical data.

Dihedral	$C_0$	$C_1$	$C_2$	$C_3$	$C_4$	$C_5$
${\phi }_1$	7.84303	0.05211	– 10.2469	0.15508	2.12905	– 0.30770
${\phi }_2$	3.51008	0.77797	– 12.7412	– 2.19178	11.5901	1.39879

Table 2. Force constants and equilibrium bond lengths and angles for the coarse-grained models used.

	$l_0/$	$k_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}/$	${\theta }_0/$	$k_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}/$
Mesogen	${\sigma }_0$	$(ϵ{\sigma }_0^{-2})$		$(ϵ\mathrm{r}\mathrm{a}{\mathrm{d}}^{-2})$
Three-site model^a	0.0	50.0	45	200
Four-site model^b	0.0	100.0	120,140	200

^aGeneric three-spherocylinder model. ^bTwo different geometeries were employed for the four-spherocylinder model.

I was in the little Italian town of Erice, perched high above the natural harbour of Trapani, on the westernmost part of Sicily, invited by Professor Claudio Zannoni to attend the International School of Liquid Crystals at the world famous Ettore Majorana centre, set in four restored monasteries. I realised that week that Claudio loved Erice (hence, the yearly International School set in such a beautiful place); but, he loved liquid crystals even more. Particularly, he loved the insights into liquid crystalline structure and dynamics provided by molecular simulation models. Over the years, the modelling work (and experiments too) from the Zannoni group has provided extraordinary insights into how liquid crystal molecules pack and arrange themselves in so many subtle ways. Consequently, I am honoured to be asked to write an article for this Festschrift and maybe it is fitting that the idea for this paper arose from discussions in Erice.

Of course, ideas born over Marsala may gloss over some of the following finer details:

• How can we produce a coarse-grained mesogenic model that is super quick, quick enough and sensitive enough to be able to provide insights into the effects we are interested in?

• How can we do sensible molecular dynamics (MD) by linking together anisotropic units without the integration of the equations of motion becoming intractable?

• How can we introduce structural chirality into the model in a way that mimics what might occur in bent-core mesogens?

• How can we measure chirality in a convenient way within a simulation model?

Many years later, all these problems are now soluble; and we are at last able to write the paper! Here, (like so many others in the field) I am again indebted to Claudio because the insights that Claudio has provided in numerous papers, over many years, have provided the knowledge and the inspiration to tackle these (and many other) problems. Grazie mille Claudio.

Mark Wilson, Durham, May 2018.

2. Introduction

The importance of chiral interactions in molecular systems has a long history dating back to the pioneering work of Pasteur. Back in 1848, Pasteur was able to separate, by hand, left-handed and right-handed crystals of tartaric acid salts to provide resolution of two enantiomeric forms [1]. Here, chiral interactions at the surface of growing crystals are sufficiently strong to cause spontaneous enantiomeric resolution and the formation of enantiopure compounds.

In liquid crystal systems, chirality is ubiquitous: from chiral mesophases to the humble twisted nematic display. Back in 1888, Friedrich Reinitzer discovered the first liquid crystal from a derivative of cholesterol [2,3]. The cholesteryl benzoate molecules that Reinitzer worked with were chiral; but, it turned out that this ‘molecular chirality’ was transmitted to the whole phase. Hence, the first liquid crystal phase discovered was itself chiral: the cholesteric (chiral nematic) phase.

Intriguingly, liquid crystals can show chirality without the molecules themselves being chiral [4]. Under the right conditions, achiral molecules can give rise to chiral structures [5–8]. This emergent chirality is a fascinating phenomenon and manifests itself in many different ways. For example, doping a cholesteric liquid crystal with an achiral bent-core molecule can lead to an increase in the helical twist (decrease in pitch) of the cholesteric helix [9,10]. Contrary to expectations, the achiral molecule, in this case, behaves like a chiral dopant with a stronger helical twisting power (HTP) than the cholesteric host. In fact, for bent-core dopant molecules, subtle changes, such as the lengthening of an achiral chain, can lead to increases in the twist of the chiral host [11,12]. We know also that achiral bent-core mesogens can form the twist-bend nematic $(N_{\mathrm{T}\mathrm{B}})$ phase [13], where molecules form a short pitch local heliconical organisation, similar to that seen in the $N_{\mathrm{T}\mathrm{B}}$ phase of achiral liquid crystal dimers [14–18] and the related $Sm{C}_{\mathrm{T}\mathrm{B}}$ phase [19]. Moreover, for achiral bent-core molecules, the presence of spontaneous saddle splay deformation is sufficient to prevent long-range orientational or translational order, leading to the formation of structures such as the dark conglomerate phase, which exhibits large homochiral domains [20–22]. Whilst spontaneous racemisation of a crystal growing from a solution of two distinct enantiomeric forms seems a beautiful natural manifestation of chiral interactions, chirality arising from strictly achiral molecules in a fluid, in which molecules are able to sample a wide range of different conformations, seems truly remarkable!

This paper seeks to provide insight into such systems from a molecular viewpoint via the use of molecular simulation models. We show that achiral bent-core molecules can possess chiral conformations that have exceptionally high HTPs. We show that the bend angle in the centre of a bent-core mesogen can influence the values of HTP and that achiral chains couple to an instantaneous chirality arising from twist within the molecular core. We demonstrate chirality transfer between molecules, leading to the preferential selection of chiral conformers of one handedness over their mirror images. We also show that nanophase segregation can occur into chiral nanodomains, mediated by preferential chiral interactions.

3. Computational

3.1. Atomistic modelling

Atomistic MD simulations were carried out using the GROMACS 4.5.5 package [23]. The calculations use our own modified version of the GAFF force field (GAFF-LCFF), which has been developed in two recent publications [24,25], including for bent-core mesogens [25].

The energy function employed in the MD simulations is given by

(1) $\begin{array}{rl}& E^{\mathrm{M}M}=\sum _{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{s}}{K}_{\mathrm{r}}(r-r_{\mathrm{e}\mathrm{q}}{)}^2+\sum _{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}\mathrm{s}}{K}_{\theta }(\theta -{\theta }_{\mathrm{e}\mathrm{q}}{)}^2\\ & +\sum _{n=0}^5{C}_n(cos(\psi ){)}^n\\ & +\sum _{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{s}}{k}_d\left(1+cos\left(n_d\omega -{\omega }_d\right)\right)\\ & +\sum _{\mathrm{i}>\mathrm{j}}^{\mathrm{N}}\left[4{ϵ}_{ij}\left({\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right)+\frac{1}{4\pi {\epsilon }_0}\frac{q_i{q}_j}{r_{ij}}\right],\end{array}$(1)

where $r_{\mathrm{e}\mathrm{q}},$ ${\theta }_{\mathrm{e}\mathrm{q}}$ are, respectively, natural bond lengths and angles, $K_r,$ $K_{\theta }$ and $C_n$ are, respectively, bond, angle and torsional force constants, ${\sigma }_{ij}$ and ${ϵ}_{ij}$ are the usual Lennard-Jones parameters and $q_i,q_j$ are partial electronic charges. Changes in $E^{\mathrm{M}\mathrm{M}}$ arising from deviations in improper dihedral angles, $\omega$, are represented by cosine functions using the force constants, $k_d$, the harmonic coefficients, $n_d$, and the phase angles ${\omega }_d$. The standard Lorentz-Berthelot mixing rules, ${ϵ}_{ij}=({ϵ}_i{ϵ}_j{)}^{1/2}$ and ${\sigma }_{ij}=({\sigma }_i+{\sigma }_j)/2$, have been applied throughout this work.

The Antechamber software from AmberTools 1.4 was used to generate GAFF topologies, with the point charges derived from the AM1-BCC method. GAFF-LCFF force field parameters from [24,25] were added in this stage. The GAFF topologies and coordinate files were converted into the GROMACS format using the ACPYPE script [26]. Atomistic simulations were carried out for one molecule in the gas phase using a stochastic dynamics (SD) methodology, a 1 fs time-step and a cutoff of 1.2 nm. Typical run lengths were in the region of 500 ns ($5\times {10}^8$ steps). Further, SD calculations employing restraints on different parts of the molecule were carried out also, as described in Section 4.1.

The original GAFF force field parameters do not accurately represent the two dihedral angles ${\phi }_1$ and ${\phi }_2$ within the benzalaniline (N-benzylideneaniline) fragment, as shown in Figure 1. Moreover, GAFF-LCFF [24,25] has not yet been extended to include these. To rectify this, we calculated torsional parameters from electronic structure theory calculations and fitted these to atomistic potentials. The quantum chemical potentials were calculated using the Gaussian 09 package [27]. The B3LYP hybrid functional was used to calculate the energies in combination with a 6–311+G* basis set. The resulting potential (after subtraction of contributions from the rest of the force field) was fitted to the Ryckaert-Bellemans function in Eq. (2) by minimisation of the squared differences between quantum chemical and force field energies. The fitted parameters are given in table 1. As shown in Figure 1, ${\phi }_1$ has minima at $0^{\circ }$ and $\pm {180}^{\circ }$ resulting in a preferred planar arrangement of the connecting group with the phenyl ring. Whilst ${\phi }_2$ prefers angles of $\pm {45}^{\circ }$ and $\pm {135}^{\circ }$ corresponding to a nonplanar arrangement. This nonplanar preference coupled with a preferred nonplanar arrangement of linking ester groups [25,28] gives rise to twisted conformations for many bent-core molecules, which, following Kawauchi et al. [28], we suggest is important in contributing to the nature of the neighbour-neighbour interactions exhibited by these molecules.

Figure 1. (Colour online) Fitted dihedrals for a benzalaniline (N-benzylideneaniline) molecular fragment.

3.2. Coarse-grained modelling

A coarse-grained simulation model was designed in which a mesogen was represented by $L/D=4$ spherocylinders of width $D={\sigma }_0=1$, joined together to form a continuous ‘bent tube’ with flexible joints, as shown in Figure 2. We choose two models: one with four spherocylinders arranged to resemble a bent-core liquid crystal molecule, and a generic model composed of three spherocylinders. The latter is chosen as the minimum number of sites needed to generate true structural chirality within a ‘rigid’ bent-core framework (i.e. a bent-molecular core linked to one additional site, which can sample different orientations relative to the bent-core, as shown in Figure 2(b)).

Figure 2. (Colour online) (a) Schematic representation of the three-spherocylinder and four-spherocylinder coarse-grained models. The diagram for the four-spherocylinder model shows the position of the four points used for the definition of the dihedral angle ${\gamma }_{1,2,2^{\prime },3}$. (b) A periodic dihedral potential around a spherocylinder is employed to favour desired equilibrium dihedral angles and provide energetic control over conformational chirality. (c) Three-dimensional representation of an instantaneously chiral conformation.

Adjacent spherocylinders were bonded together through harmonic springs, $U_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}=k_{\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}}(l-l_0{)}^2$, linking the ends of the spherocylinder line segments. Additional harmonic potentials, $U_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}=k_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}(\theta -{\theta }_0{)}^2$, were used to constrain the angles of terminal spherocylinders with respect to the axis vector of the central spherocylinder for the three-site model and with respect to neighbouring spherocylinders for the four-site model. Equilibrium values for bond lengths and angles are given in Table 2. For this phenomenological model, which is akin to a DPD model for a liquid crystal mesogen [29] (but with anisotropic sites), force constants are chosen to be sufficiently large to maintain the molecular geometry for the chosen simulation conditions and the soft-core spherocylinder potential described in the following, whilst sufficiently small to conserve energy for the chosen time step, $\delta t$. As usual, nonbonded interactions were excluded between directly bonded spherocylinder sites.

To control the conformational chirality of the model molecule, a ‘dihedral’ potential was employed to control rotation around a spherocylinder, as defined by the angle ${\gamma }_{ijkl}$ (Figure 2(b)). Two dihedral potentials were considered,

(2) $U_{\mathrm{d}\mathrm{i}\mathrm{h}}=A\left[1+cos\left(n{\gamma }_{ijkl}-{\gamma }_0\right)\right]$(2)

and

(3) $U_{\mathrm{d}\mathrm{i}\mathrm{h}}=-A\left[cos{\left(\frac{{\gamma }_{ijkl}}{2}-{\gamma }_0\right)}^n+cos{\left(\frac{{\gamma }_{ijkl}}{2}+{\gamma }_0\right)}^n\right]+A$(3)

Here, the parameter $A$ controls the height of the barrier between the consecutive wells and ${\gamma }_0$ is used to change the equilibrium dihedral angle (and hence the chirality). The integer $n$ sets the number of wells in Eq. (2) and controls the width of the two wells in Eq. (3). In our calculations, for ${\gamma }_{ijkl}$ we measure the torsional angle between ends of spherocylinder line segments for a central spherocylinder and the centres of mass of the two adjacent spherocylinders (i.e. ${\gamma }_{ijkl}={\gamma }_{1,2,2^{\prime },3}$ would be the single dihedral defined for the three-spherocylinder model) and integrate the forces and torques that arise from the chosen potential (Eqs. (2) or (3)). For the three-site model, there is just one dihedral angle defined, ${\gamma }_{ijkl}={\gamma }_{1,2,2^{\prime },3}$. However, for the four-site model, there are two dihedral angles: ${\gamma }_{1,2,2^{\prime },3}$ and ${\gamma }_{2,3,3^{\prime },4}$ as shown in Figure 2. Note that, for the three-spherocylinder model, choosing values of ${\gamma }_{ijkl}={\gamma }_{1,2,2^{\prime },3}={180}^{\circ }$ or ${\gamma }_{1,2,2^{\prime },3}=0^{\circ }$ results in an achiral (planar) molecule; and, for the four-spherocylinder model, choosing both ${\gamma }_{1,2,2^{\prime },3}$ and ${\gamma }_{2,3,3^{\prime },4}$ equal to values of either ${180}^{\circ }$ or $0^{\circ }$ similarly results in an achiral (planar) molecule.

Nonbonded interactions between spherocylinders were described by a soft-core spherocylinder model [30,31] (see in the following) with a parametrisation, $U_{\mathrm{m}\mathrm{a}\mathrm{x}}=100.0{ϵ}_0,U_{\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{r}}=2000.0{ϵ}_0,{ϵ}_1=220.0{ϵ}_0,$ ${ϵ}_2=0.0{ϵ}_0$. This parameterisation is known to lead to nematic stability over a wide density range for single site mesogens [30,31]. We employed the MD methodology of reference [30] using the GBMOL program [32–34]. Simulations were carried out at $T^{\ast }=kT/{ϵ}_0=3$ in the nematic phase, with typical runs employing 2–4 million MD steps.

The soft-core spherocylinder interaction potential is given by

(4) $U^{\ast }(d)=\left\{\begin{array}{cc}{U}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\ast }{(1-d^{\ast })}^2+{ϵ}^{\ast },& \mathrm{i}\mathrm{f}{d}^{\ast }<1\\ U_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\ast }{(1-d^{\ast })}^2-U_{\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{r}}^{\ast }({\hat{r}}_{ij},{\hat{e}}_i,{\hat{e}}_j){(1-d^{\ast })}^4+{ϵ}^{\ast },& \mathrm{i}\mathrm{f}1\le d^{\ast }<d_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\ast }\\ 0,& \mathrm{i}\mathrm{f}{d}^{\ast }\ge d_{\mathrm{c}\mathrm{u}\mathrm{t}}^{\ast }\end{array}\right.$(4)

where $U_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\ast }$ is the interaction energy of fully overlapped particles, $d^{\ast }$ is the line segment distance and ${ϵ}^{\ast }$ is the well depth. $U_{\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{r}}^{\ast }({\hat{r}}_{ij},{\hat{e}}_i,{\hat{e}}_j)$, the magnitude of the attractive component is angle dependent, i.e. a function of the orientation of the vector from a centre of mass, ${\hat{r}}_{ij}$, and orientations of particles i and j (${\hat{e}}_i$ and ${\hat{e}}_j$ ). This potential can be easily tuned for different liquid crystal behaviours by altering the four parameters: $U_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\ast }$, $U_{\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{r}}^{\ast }$, ${ϵ}_1$ and ${ϵ}_2$ [30,31] in a similar way to the well-known Gay-Berne potential [35].

For the three-site spherocylinder model, $N=4000$ molecules were simulated in a cuboidal box ($L_x=L_y\approx 41.35{\sigma }_0$ and $L_z\approx 58.48{\sigma }_0$) with periodic boundaries in $x$ and $y$ directions. A repulsive wall was employed in the $z$ -direction, acting on the $z$ -component of the centre of mass of spherocylinders, $r_z$, beyond the top or bottom of the simulation box: $U_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}=A_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}(1-r_z/{\sigma }_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}{)}^2$ with $A_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}=100.0{ϵ}_0$ and ${\sigma }_{\mathrm{w}\mathrm{a}\mathrm{l}\mathrm{l}}=3{\sigma }_0$.

3.3. Importance of choosing soft-core interaction potentials

The key advantage of the coarse-grained DPD-like model used in this work is the speed of equilibration. The model is sufficiently repulsive to promote liquid crystal behaviour through anisotropic excluded volume interactions but sufficiently soft to allow rapid movement through phase space, avoiding the jamming which is sometimes seen with harder particle coarse-grained models. The benefits of this type of model have been demonstrated by Zannoni and co-workers [36,37] and used to look at changes in molecular organisation in the actuation cycle of a liquid crystal elastomer [38]. We also note the similarity of our four-site model to the Monte Carlo model of Peroukidis et al. [39]. In their model, they used spherocylinders of different sizes for the core and ‘tail’ of a bent-core mesogen, allowing the molecule to adopt both chiral and achiral conformations. They observed the formation of clusters exhibiting local smectic order, which were either orthogonal or tilted, together with homochiral domains through a helical arrangement of the tilted smectic clusters.

We note in passing that the use of a soft repulsive wall together with our soft-core spherocylinder model leads to very weak anchoring at the wall surface by design. This means that the wall does not restrict the rotation of molecules in the $x-y$ plane, allowing a twist induced on a nematic director in the $x-y$ plane to easily overcome any influence that can otherwise be introduced by the presence of harder walls. In addition, the weak anchoring ensures that smectic layers are not induced at the wall surface. An interesting alternative will be the use of soft walls in combination with self-determined boundary conditions, as in the recent work of Kuhnhold and Schilling [40].

3.4. Quantifying molecular twist in the bulk

To quantify the effect of molecular chirality on the twist of the bulk phase, for boxes with walls in the $z$ -direction, the simulation box was divided into layers with thickness of 1.7${\sigma }_0$ along the $z$-axis. A local nematic director was calculated for each layer, $n_{\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}$ from diagonalisation of the ordering tensor for the central spherocylinders. This allowed the calculation of the twist angle, $\varphi (z/{\sigma }_0)$, along the $z$-axis. For linear increases in $\varphi$ along $z$, the bulk pitch, $p$, can be obtained by a linear fit, $\varphi (z^{\ast })=a{z}^{\ast }+b$, as $p=2\pi /a$.

3.5. Calculation of HTPs

Ferrarini et al. [41–44] showed that molecular isosurfaces could be used to evaluate the chirality order parameter, $\chi$, which is proportional to the HTP of a mesogen in a nematic phase. To calculate $\chi$, we use the step-by-step approach described in detail by Earl and Wilson [45]. Specifically, we first calculate a molecular surface from an atomistic molecular conformation (or a coarse-grained representation in the case of the coarse-grained models). We used an isosurface corresponding to the van der Waals surface accessible to a spherical probe with a diameter of 5 Å (in the case of coarse-grained models the probe was placed in contact with the surface of the spherocylinders). Isosurfaces were generated using the simple invariant molecular surface routine of Vorobjev and Hermans [46] with a resolution of 10 dots per Å${ }^2$. We next calculated a normal vector at each surface point and then calculated the surface tensor $T$, the helicity tensor $Q$ and the ordering matrix $S$ (see Eqs. (3), (4) and (6) in [45]) in order to evaluate $\chi$. $\chi$ values were determined for $1\times {10}^6$ conformations for atomistic simulations and for selected snapshots for coarse-grained models.

4. Results and discussion

4.1. Conformational distribution functions

Angular distribution functions, obtained from the SD simulations, are plotted in Figure 3 for the bend angle of the molecular core and the angle between the core and chain. As shown in Figure 3, the measured distributions are highly dependent on the exact definition of the angle. However, in the most commonly used definition of the ‘molecular core’ ($A_1$ in Figure 3) the angle distribution function is peaked at 120${ }^{\circ }$ commensurate with that expected from a bent structure originating from a central phenyl ring. The distribution function for $A_1$ is, however, quite broad, reflecting the flexibility of the ester and benzalaniline linkages in the core. Consequently, the idealised concept of a ‘rigid bent-core’ is certainly not true. The core-chain angle has a particularly wide distribution function, which is peaked between ${140}^{\circ }-{150}^{\circ }$ but which extends to quite low angles arising from a long-flexible chain, which is able to bend back on itself. Within a nematic phase, we would expect these distribution functions to be somewhat less broad, with the angle distributions being coupled to the field of surrounding molecules. The phenomenon is known well for chain conformations in $n$-alkylcyanobiphenyls, where preferential selection occurs for conformations that enhance alignment with the nematic mean field [47,48].

Figure 3. (Colour online) Angle distributions for P8OPIMB from single-molecule stochastic dynamics calculations at 300 K using the GAFF-LCFF force field.

4.2 χ values for PnOPIMB

Recently, Kim et al. [12] have calculated $\chi$ values for various members of the PnOPIMB series. Here (Figure 4), we repeat that work for $n=8,12,16$ using the new GAFF-LCFF force field, and extend the calculation to the core only $n=0$ homologue. In addition, we carry out SD calculations with MD restraints applied to the chain (giving an all-trans chain), to the ester group (giving a planar arrangement for the three central rings) and to the whole core (giving a planar arrangement for the whole core). Several features are apparent from Figure 4:

Figure 4. (Colour online) Top: $\chi$ distribution functions (arbitrary units) for a PnOPIMB molecule from stochastic dynamics calculations at 300 K. Middle: $\chi$ distribution functions for P12OPIMB employing chain, ester and core restraints. Bottom: An instantaneous twisted conformation for P8OPIMB with $\chi =806.24$ showing the molecular surface used to calculate $\chi$.

• Excellent conformational sampling within the SD calculations leads to distribution functions that are symmetrical and centred on zero; demonstrating, as expected, no net chirality.

• The wings of the distribution function extend to exceptionally high values of $\chi$, around 50$\times$ larger than values of $\chi$ that are associated with very high HTP chiral dopants, such as bridged biaryls and helicenes derivatives [45].

• The effect of the chain length is significant, with increased chain lengths leading to a widening of the distribution function. In contrast, the calculations with no chain give a sharply peaked distribution. (We note that with GAFF-LCFF, which has an improved description of the ester and benzalaniline linkages, and the alkoxy chain, that a widening of the $\chi$ distribution is seen for $n>12$, which was not observed by Kim et al.)

• For P12OPIMB, applying restraints to the whole core (giving an approximately planar core), leads to a big reduction in the higher values of $\chi$. However, applying restraints to the ester linkage and the chain leads to only small changes in the $\chi$ distribution.

We suggest that the huge instantaneous values of $\chi$ seen in the wings of the distribution are responsible for the increase in HTP when a chiral nematic is doped with PnOPIMB. Highly twisting conformations will be preferentially selected in a chiral molecular field. Moreover, it appears that the main origin of the high HTP values arises from twisted (chiral) conformations of the core and if these are eliminated, by using a planar core, the $\chi$ distribution is sharpened up considerably. We note, however, that even in this case, chiral conformations of the chains can occur and their HTP is amplified by the presence of the bent-core shape. Constraining the chains to an all-trans conformation does not reduce HTP significantly for P12OPIMB. This is because the twist of the core still leads to an overall twisted shape for the molecule, with the two chains pointing in different directions, and hence a large HTP is obtained. Interestingly, if the chains are extended to $n=16$, applying chain restraints leads to a further widening of the distribution function. This suggests that, in this case, it is not the addition of further ‘instantaneously chiral’ chain conformations that is most important but rather the length of the chain that contributes to an overall enhanced chiral shape for the molecule as a whole. We explore these ideas further with the use of idealised coarse-grained models (see the following).

4.3 χ values for coarse-grained models

We carried out single-molecule dynamics calculations for the three-site and four-site models, described in Section 3.2, using an equilibrium bond angle of ${120}^{\circ }$ and an equilibrium dihedral of ${180}^{\circ }$ but with a very low force constant ($k=10\epsilon$) to allow easy sampling of twisted conformations arising from other values of this dihedral. In addition, we removed the two terminal spherocylinders from the four-site model, leaving only a ‘bent-core’ ($n=2$ model). The results are given in Figure 5(a).

Intriguingly, the $n=3$ system appears to be the model that exhibits chiral conformations with the largest HTP. However, the four-site model captures the same shaped $\chi$ distribution that we saw from atomistic modelling. Not surprisingly, removing the two terminal spherocylinders to leave a $n=2$ system leads to a distribution strongly peaked about $\chi =0$ with only configurations where the spherocylinders fluctuate out of plane contributing significantly to $\chi$. For the $n=3$ and $n=4$ models, the $\chi$ values in the wings of the distribution are an order of magnitude larger than those obtained with the highest HTP conformations seen in the atomistic simulations. This provides hope that, for these models, one might be able to see induced helical twists on the scale of molecular simulations (see Section 4.4).

Further insights are provided by examining how $\chi$ varies with molecular conformation. For the three-site model (Figure 5(b)) the value of $\chi$ peaks for a dihedral of ${\gamma }_{1,2,2^{\prime },3}={105}^{\circ }$ $(\theta ={120}^{\circ })$ and varies smoothly as a function of this dihedral. The values of $\chi$ are also shown to depend on the internal angle $\theta$. For the four-site model, with internal angles of 120${ }^{\circ }$, the highest HTP occurs at ${\phi }_1=-{\phi }_2={90}^{\circ },-{90}^{\circ }$, and there is an unexpected dip in HTP for ${\phi }_1=-{\phi }_2={45}^{\circ },-{45}^{\circ }$ (Figure 5(c)). Moreover, changing the internal angle to ${140}^{\circ }$ (Figure 5(d)), as would occur in an oxadiazole-based bent-core mesogen [49], subtly changes the form of the curve, with the maxima shifting to ${\phi }_1=-{\phi }_2={75}^{\circ },-{75}^{\circ }$ and the local minima shifting to ${\phi }_1=-{\phi }_2={30}^{\circ },-{30}^{\circ }$. We note that, for the four-site model, conformations where ${\phi }_1={\phi }_2$ (which, for an achiral molecule, always have identical torsional energies to conformations with ${\phi }_1=-{\phi }_2$) have superimposable mirror images and hence $\chi =0$.

Figure 5. (Colour online) $\chi$ values for single molecule coarse-grained models.

4.4 Bulk simulations of the three-site coarse-grained model

The effect of molecular-level chirality on a bulk phase was initially determined via six simulations using the dihedral potential of Eq. (2); using parameters $A=20.0{\epsilon }_0,n=1$ and ${\gamma }_0=0^{\circ }$, ${10}^{\circ }$, ${30}^{\circ }$, ${45}^{\circ }$, ${60}^{\circ }$, ${70}^{\circ }$. These correspond to equilibrium dihedrals ${\gamma }_{\mathrm{e}q}={\gamma }_{1,2,2^{\prime },3}={180}^{\circ }$, $-{170}^{\circ }$, $-{150}^{\circ }$, $-{135}^{\circ }$, $-{120}^{\circ }$, $-{110}^{\circ }$. The simulations used the soft repulsive wall boundary conditions described in Section 3.2, leading to macroscopic twists across the simulation box, corresponding to a chiral nematic phase with a small pitch length. We quantify the twist by calculating the order parameter, $⟨S_2⟩=⟨P_2(cos\theta ){⟩}_{\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}$, and the director, $n_{\mathrm{l}ayer}$, along the $z$ -axis of the simulation box (as described in Section 3.4), and by plotting the change in angle of the director with changing $z$, as shown in Figure 6. The overall twist, $\varphi$, is seen to increase from, $\varphi \approx 0^{\circ }$, for the planar (achiral) molecules $({\gamma }_{\mathrm{e}\mathrm{q}}={180}^{\circ })$ to $\varphi \approx {278}^{\circ }$ for the highest chiral strength $({\gamma }_{\mathrm{e}\mathrm{q}}=-{110}^{\circ })$. The mean nematic order parameter remains well inside the nematic region, ranging from $⟨S_2⟩\approx 0.67\pm 0.02$ for the achiral system to $S_2\approx 0.604\pm 0.004$ for the system with the highest chirality. A small local increase in $S_2$ near the wall was observed (maximum of 17%) which drops to the bulk value over the width of one histogram bin $(1.7{\sigma }_0)$.

Figure 6. (Colour online) Top: Snapshots for $N=4000$ for the three-site spherocylinder model. Snapshots are viewed along the $z$ -axis of the simulation box for two chiral nematic systems, with ${\gamma }_{\mathrm{e}q}=1500$ (top left) and ${\gamma }_{\mathrm{e}q}=1200$ (top right). Bottom: The effect of molecular chirality on the twist angle between a reference nematic director, $n_{\mathrm{r}\mathrm{e}\mathrm{f}}$, at the base of the simulation box and the local nematic director, $n_{\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}}$, along the $z$-axis of the simulation box.

To model the effects of adding achiral dopants (with the possibility of chiral conformations) to a system, achiral molecules were designed with two symmetrical potential wells at ${\gamma }_{\mathrm{e}\mathrm{q}}=\pm {110}^{\circ }$, separated by a barrier of height $20{\epsilon }_0$ (corresponding to $A=20\epsilon ,n=40$ and ${\gamma }_0={55}^{\circ }$ in Equation (3)). Two host systems with 10% doping were considered: a chiral host phase with $({\gamma }_{\mathrm{e}\mathrm{q}}=-{150}^{\circ })$ and, as a control (to check statistics) a second host achiral host phase $({\gamma }_{eq}={180}^{\circ })$ was also employed. Both cases used a maximum barrier height of $40{ϵ}_0$.

Figure 7 shows the dihedral angle distributions, $S({\gamma }_{ijkl})$, obtained for the solvent and dopant molecules. The preferred handedness of the dopant molecules can be obtained from the relative populations of right- and left-handed molecules, $n_L$ and $n_R$, provided by numerical integration of $S({\gamma }_{ijkl})$. For the test case, i.e. dopant molecules in the achiral nematic phase (dotted line in Figure 7(a)), $n_L\approx n_R$ within a statistical error of 1%. This level of sampling is only possible because of the speed of the coarse-grained model we employ, which means that it is able to fully equilibrate the system and sample the conformations to a high degree of statistical accuracy.

When the same dopants are added to the chiral host phase (dot-dashed line in Figure 7(a)) preferential selection of conformations occurs for the conformer with the same handedness as the host phase. Here, relative populations of $n_{\mathrm{L}}\approx 59.3\pm$ 0.5% and $n_{\mathrm{R}}\approx 40.7\pm 0.5$ were observed for a left-handed host solvent. This same observation is found for a series of values of different ${\gamma }_{\mathrm{e}q}$ for the dopant system. In the case of ${\gamma }_{\mathrm{e}\mathrm{q}}=-{150}^{\circ }$ , the dopants have conformations with higher twisting power than the host phase. A subsequent plot of bulk twist angle in Figure 7(b) shows an increase in the twist of the system with 10% achiral dopants in comparison to the pure chiral solvent. From the linear fits, pitch lengths were estimated as $p=(88.5\pm 1.0){\sigma }_0$ and $p=(103.2\pm 0.7){\sigma }_0$ for the doped and pure systems, respectively. This corresponds to a 17% decrease of the pitch length (increase in the twist).

Figure 7. (Colour online) (a) Dihedral angle distributions from MD simulations. Chiral solvent (solid line), dopant in chiral solvent (dashed line), dopant in nematic solvent (dotted line) and solvent without dopant molecules (dot-dashed line). (b) Twist angle for nematic director along the $z$ -axis: pure melt (squares) and system with 10% doping of achiral solute. The lines are least-squares fits of the linear function $f(z)=az+b$. (c) Radial distribution functions, $g(r)$, for solvent molecules (solid line), dopant molecules (dashed line) and solvent-dopant cross distribution (dotted line). (d) Snapshot of the dopant molecules, where solvent molecules are removed for clarity.

Three possible mechanisms can be evoked to explain the increase in the bulk twist of the doped sample.

(1) Dopant molecules could increase the preferred dihedral angles of the solvent molecules, thereby increasing chirality of the solvent molecules.

(2) Dopant molecules are distributed evenly in the system but the preferential selection of high twisting conformations of the same handedness as the solvent leads to an overall increase of the bulk twist.

(3) The dopant molecules could undergo microphase separation to form highly twisting clusters leading to an overall increase in the twist across the system.

The first of these mechanisms can be immediately ruled out. The calculated dihedral angle distributions for the solvent with 10% doping (solid line in Figure 7) and the pure solvent (dot-dashed line) are in agreement within the line width. We expect also that this mechanism is not experimentally important, i.e. most chiral hosts are chiral due to the presence of chiral centres or in-built structural chirality, and do not tend to have highly twisting conformations, which would come into play only when preferentially selected through interaction with an achiral dopant. From the results in Figure 7, the preferential selection of high HTP conformers of the dopant mechanism clearly does occur. However, this might be seen in conjunction with the third mechanism. To test this, we looked at the spatial distribution of dopant molecules by calculating partial radial distribution functions, $g(r)$ for all pairs of species (Figure 7(c)). Each $g(r)$ shows a peak at a separation of $r\approx 1.3{\sigma }_0$ corresponding to nearest neighbour interactions. However, the $g(r)$ for the dopant-dopant interactions shows the growth of a slight shoulder between separations of $r\approx 2.4{\sigma }_0$ and $r\approx 6{\sigma }_0$, indicating a slight preference for the formation of small clusters at this distance. A snapshot picture (Figure 7(d)) confirms that at the 10% dopant level the dopant molecules are uniformly distributed across the simulation box, apart from some evidence for the formation of very small local clusters. Hence, it appears that preferential interaction between chiral conformations of neighbouring dopants is likely to contribute to an increase in overall HTP.

Further evidence of chiral selectivity can be observed in related simulations. In Figure 8, we show a snapshot picture taken from a uniform nematic phase for the three-site spherocylinder model. Here, a system of strictly achiral molecules (with a preference for no chirality ensured from ${\gamma }_{\mathrm{e}q}={180}^{\circ }$) has been simulated and a colour coding has been applied according to the handedness of the molecules. From this, it can be seen how locally the molecules prefer the same handedness as their neighbours leading to nanophase segregation into small domains with opposite handedness. In our model, within a nematic phase, the chiral domains are quite small as opposed to experimental findings in more complex phases [20,21,50,51] where large homochiral domains are sometimes found, enhanced by strong neighbour packing effects. Nonetheless, the preference for chiral interactions seen in our model is clear.

Figure 8. (Colour online) A snapshots showing chiral nanophase segregation in an achiral nematic phase (${\gamma }_{\mathrm{e}\mathrm{q}}={180}^{\circ }$) where the colour coding has been applied as blue (red) left (right) handed conformations.

5. Conclusions

We have studied chirality arising from achiral bent-core liquid crystal molecules via atomistic and coarse grained-simulation models. We demonstrate huge HTPs for individual conformations of three bent-core molecules, PnOPIMB ($n=8,12,16$), with the HTP values increasing with alkoxy chain length. We show that these conformational effects can be captured using simple three-site and four-site soft spherocylinder models. Further, we show, with a three-site spherocylinder model, that chiral induction can take place within a bulk nematic phase.

For achiral molecules our coarse-grained simulations show a uniform nematic phase; but, when conformational chirality is introduced, a chiral nematic phase is formed with the measured twist angle of the bulk phase increasing as the internal twist angle of individual molecules is increased. Moreover, the simulations demonstrate a decrease in pitch length (increased twist) for a chiral nematic, when the system is doped with an achiral dopant, in agreement with surprising experimental results for bent-core molecules [9,11]. Here, we see the preferential selection of high twisting enantiomeric conformations over their mirror images.

Finally, we note that we see evidence for local chiral segregation into nanodomains of opposite handedness within a bulk nematic phase. It will be interesting to explore such effects further with a view to understanding the formation of helical nanofilaments and larger chiral domains within more exotic liquid crystal mesophases.

Acknowledgments

MRW and MW would like to thank EPSRC for funding under grants EP/J004413/1 and EP/P007864/1. JSL would like to acknowledge funding from an IdEx (Initiative d’Excellence) Bordeaux Junior Chair. We would like to thank Professor Claudio Zannoni for valuable discussions over many years, and for his insights into the simulation of liquid crystal systems. We gratefully acknowledge conversations with Professors Vanakaras and Photinos in Erice and elsewhere, which have provided many insights into the behaviour of bent-core systems.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Engineering and Physical Sciences Research Council [EP/J004413/1,EP/P007864/1].