Molecular Simulation Approaches to the Study of Thermotropic and Lyotropic Liquid Crystals

Abstract

Over the last decade, the availability of computer time, together with new algorithms capable of exploiting parallel computer architectures, has opened up many possibilities in molecularly modelling liquid crystalline systems. This perspective article points to recent progress in modelling both thermotropic and lyotropic systems. For thermotropic nematics, the advent of improved molecular force fields can provide predictions for nematic clearing temperatures within a 10 K range. Such studies also provide valuable insights into the structure of more complex phases, where molecular organisation may be challenging to probe experimentally. Developments in coarse-grained models for thermotropics are discussed in the context of understanding the complex interplay of molecular packing, microphase separation and local interactions, and in developing methods for the calculation of material properties for thermotropics. We discuss progress towards the calculation of elastic constants, rotational viscosity coefficients, flexoelectric coefficients and helical twisting powers. The article also covers developments in modelling micelles, conventional lyotropic phases, lyotropic phase diagrams, and chromonic liquid crystals. For the latter, atomistic simulations have been particularly productive in clarifying the nature of the self-assembled aggregates in dilute solution. The development of effective coarse-grained models for chromonics is discussed in detail, including models that have demonstrated the formation of the chromonic N and M phases.

1. Introduction

In recent years, molecular simulation has become a powerful tool for studying a range of soft matter systems. Simulations have been extended to bulk polymers [1], polymer surfaces [2], proteins [3,4], membranes [5,6,7], self-assembly in solution [8], and molecules at water interfaces [9,10], in addition to many other systems. In these cases, simulations aim to provide quantitative predictions to compare with experiment and also to provide qualitative insights into local molecular structure and order, which are often difficult to obtain by experimental means. These comments are particularly true for liquid crystalline systems, where an important part of the contemporary picture of how molecules are ordered within liquid crystal phases comes from the insights that molecular simulation has been able to offer over the last few decades [11,12,13,14].

Molecular simulation models are part of a traditional hierarchy of simulation methodologies which cover different time and length scales from the microscopic to mesoscopic to continuum scales. For liquid crystals, this hierarchy is particularly significant (see Figure 1) because it covers the following:

• The quantum mechanical regime, where single-molecule calculations are valuable in determining molecular properties of single thermotropic mesogens to use, for example, in the development of materials for displays [15];
• The atomistic regime where atomistic simulation can, in principle, be used to study the detailed molecular structure of a liquid crystal and predict bulk properties for simpler thermotropic phases, such as nematics [16,17];
• The coarse-grained regime where simulations can be used to study the structure of smectics [18], twist grain boundary phases [19], and polymer liquid crystals [20,21,22,23] and more complex liquid crystal phases where molecular shape and packing (or often the shape of a larger object such as a colloid) are often significant in determining the structure of the phase;
• The continuum regime where simulations no longer consider a molecular description but instead considers a fluid description where the local orientation of the director can be followed in complex geometries [24].

Over the last decade, the possibilities for modelling liquid crystals have been extended hugely with the realisation of readily accessible parallel computer architectures and simulation algorithms that can take advantage of them. This has led to a wide range of new studies [11]. This perspective article, which is part of a collection of articles showcasing liquid crystal research in the UK, highlights some of the recent progress made in studying thermotropic and lyotropic liquid crystal systems using molecular simulation methods. It covers both atomistic simulations and simulations carried out using different types of coarse-grained molecular simulation models. This article concentrates on new developments and also points to areas where the continued increase in computer power is opening up new possibilities.

2. Thermotropic Liquid Crystals

2.1. Atomistic Simulations of Thermotropics: Towards the Prediction of Accurate Transition Temperatures

A holy grail of atomistic simulation is to predict the phase sequence, transition temperatures, and local molecular structure of liquid crystal phases. However, this goal is extraordinarily difficult to realise and, in many ways, is one of the most challenging goals in soft matter simulation. Thermotropic liquid crystalline systems are often sensitive to small details in chemical structure. The addition of an extra functional group, or even just a small change in the length of an alkyl chain, can dramatically change transition temperatures or even change the sequence of thermodynamically stable liquid crystal phases that a molecule exhibits. Simulations have helped demonstrate the reasons for this but has not fully solved the prediction problem!

Firstly, the stability of thermotropic liquid crystals is a mix of two factors: shape and attractive interactions. For example, a small increase in molecular length will increase excluded volume effects and will drive molecular alignment of a thermotropic mesogen in a way that is better-known for colloidal systems and is well-understood by Onsager theory [25,26,27]. Moreover, an increase in density of packing (strictly a decrease in free volume) for a thermotropic system drives molecular alignment by a similar mechanism. However, changes in functional groups can disrupt molecular packing and also add preferred local interactions that can alter the phase behaviour. Add this to the further complication that many liquid crystal molecules can potentially exhibit a range of different unseen phases (i.e., unseen under standard thermodynamic conditions because they are marginally higher in free energy than the phases which are seen) and the phase diagram prediction problem becomes exceedingly challenging.

Nonetheless, considerable progress has been made in producing high-quality atomistic models that can represent structurally simpler liquid crystal phases and predict transition temperatures within 5–10 ${}^{°}$C [28]. The key to the successes has been improvements in the accuracy of molecular force fields.

A typical molecular force field takes the following functional form:

$$\begin{array}{ccc}\hfill E^{\mathrm{MM}}& =& {\displaystyle \sum _{\mathrm{bonds}}}{K}_{\mathrm{r}}{(r-r_{\mathrm{eq}})}^2+{\displaystyle \sum _{\mathrm{angles}}}{K}_{\theta }{(\theta -{\theta }_{\mathrm{eq}})}^2\hfill \\ \hfill & +& {\displaystyle \sum _{\mathrm{torsions}}}{\displaystyle \sum _{n=0}^5}{C}_n{(cos\left(\psi \right))}^n+{\displaystyle \sum _{\mathrm{impropers}}}{k}_d\left(1+cos\left(n_d\omega -{\omega }_d\right)\right)\hfill \\ \hfill & +& {\displaystyle \sum _{i>j}^N}\left[4{\epsilon }_{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],\hfill \end{array}$$

(1)

where $r_{\mathrm{eq}}$ and ${\theta }_{\mathrm{eq}}$ are, respectively, natural bond lengths and angles, $K_{\mathrm{r}},$ $K_{\theta }$, and $C_n$ are, respectively, bond, angle, and torsional force constants, $\psi$ is a dihedral (torsional) angle, $k_d$ is a force constant associated with improper dihedral angles, $\omega$, $n_d$, and ${\omega }_d$, respectively, represent an improper dihedral angle, its periodicity, and the phase angle associated with it. ${\sigma }_{ij}$ and ${\epsilon }_{ij}$ are the usual Lennard–Jones parameters modelling nonbonded interactions, and $q_i,q_j$ are partial electronic charges. Here, we have written $E^{\mathrm{MM}}$ in the form used by the popular General AMBER Force Field (GAFF) [29].

In terms of liquid crystal systems, several things are crucial:

• A good charge distribution to represent the local electrostatic potential around a thermotropic mesogen;
• Accurate torsional potentials to represent internal rotations, as these can dramatically alter the average shape exhibited by a thermotropic mesogen within a liquid crystalline phase;
• Excellent nonbonded interactions to represent steric repulsion and attractive interactions.

In fact, the torsional potentials and nonbonded interactions are absolutely critical in determining transition temperatures. The majority of mesogens have alkyl chains that are “melted” in a mesophase (i.e., not in the lowest energy conformation but exhibit a range of conformations). If the torsional interactions are, for example, too stiff, then the average structure will tend to be elongated, phases will tend to be artificially stabilised; as a result, transition temperatures will be increased. Likewise, a very small error leading to an increase in density can dramatically promote mesophase stability by excluded volume effects and the Onsager mechanism [30], i.e., promoting translational entropy at the expense of rotational entropy.

With these factors in mind, Boyd and Wilson developed a version of GAFF, GAFF-LCFF (GAFF-Liquid Crystal Force Field), that was tuned for liquid crystalline systems [28,31]. Their work involved the careful optimisation of torsional potentials by fitting to high-quality density functional torsional scans, and (following the approach of previous workers [32,33,34]), careful optimisation of nonbonded parameters to reproduce densities and heats of vaporization of small-molecule fragments. Here, in particular, they aimed, where possible, to fit the densities of small molecules to better than 1%. Whilst such work has traditionally been extremely time consuming, in the future, it is likely to become much easier through developments such as the Open Force Fields project and automatic fitting tools such as the ForceBalance software [35].

Using GAFF-LCFF, liquid crystal state points and the isotropic to nematic transition are accessible for simulations extended over ∼240 ns per state point (although simulations times required are strongly system size-dependent and also depend on the viscosity of the system). These time scales are very short in comparison to director reorientation in bulk liquid crystals but on the nanoscale are sufficient to see the spontaneous alignment of a nematic liquid crystal from the isotropic phase. Figure 2 shows snapshots illustrating the molecular structure in the nematic phase of the bent-core mesogen C5-Ph-ODBP-Ph-OC12, including the beginnings of microphase phase separation between core and chains.

Orientational order parameters are typically calculated from suitable vectors within the molecule or, in the case of many calamitic mesogens, from the molecular long axis obtained from the moment of inertia tensor. The instantaneous average across molecular vectors in the simulations leads to both a liquid crystal director $\overrightarrow{n}$ and an orientational order parameter $S_2$. In practice, these are obtained by calculating the ordering tensor:

$${\mathbf{Q}}_{\alpha \beta }\left(t\right)=\frac{1}{2N}\sum _{i=1}^N\left[3u_{i\alpha }{u}_{i\beta }-{\delta }_{\alpha \beta }\right],\alpha ,\beta =x,y,z,$$

(2)

where the sum runs over all N molecules. The largest eigenvalue of the $\mathbf{Q}$ tensor represents the following:

$$P_2\left(t\right)=\frac{1}{N}\sum _{i=1}^N{P}_2(cos{\theta }_i),$$

(3)

where $P_2$ is the second Legendre polynomial, and the associated eigenvector is the director $\overrightarrow{n}\left(t\right).$ $S_2$ is defined as the time average of $P_2$ over a suitable time interval where $〈P_2〉$ is unchanging. However, in practice, to minimise system size effects in locating the phase transitions, $S_2$ is often obtained from $-2\times$ and the middle eigenvalue of $\mathbf{Q}$, which fluctuates about a value of zero in the isotropic phase but equals $P_2\left(t\right)$ in the nematic phase. Here, we note in passing that in relatively small simulated systems, there is always a danger of some uncertainty connected with system size and the possibility of supercooling past phase transitions. Although this is less of a problem with nematic liquid crystals than many other soft matter systems where the enthalpy change associated with the phase transition is larger.

For the systems tested (molecular structures given in Figure 3), GAFF-LCFF provides excellent liquid crystal clearing point predictions, as shown in

Table 1. Experimental and simulated clearing points for a series of liquid crystalline systems using GAFF-LCFF modifications [28,31] relative to the GAFF force field [29].

Molecule	${\mathit{T}}_{\mathbf{NI}}$ (exp.) /K	${\mathit{T}}_{\mathbf{NI}}$ (GAFF-LCFF) /K
1,3-benzenedicarboxylic acid,1,3-bis(4-butylphenyl)-ester	452	450–460
C5-Ph-ODBP-Ph-OC12	512.6	∼510
C4-Ph-ODBP-Ph-C7	507	∼500
C4O-Ph-ODBP	558	550–560
C4O-Ph-ODBP (trimethylated)	421	420–430

. Because of its success, GAFF (with GAFF-LCFF modifications) has been used for a range of liquid crystal systems, including recent studies into de Vries behaviour [36,37].

2.2. Simulation Insights into the Structure of New Phases

One of the most valuable features of atomistic simulations is that they can provide some insight into the formation of new liquid-crystalline phases and shed some light on liquid–crystalline structures where the arrangement of molecules is unclear from experiments. Figure 4 is from the preliminary work of Yu and Wilson and shows the structure of the twist-bend nematic (${\mathrm{N}}_{\mathrm{TB}}$) phase for liquid crystal dimer CB7CB. Here, CB7CB is gradually cooled from a nematic phase into the underlying phase and equilibrated over periods of 100 ns per state point. There was originally considerable debate about the structure of the phases of CB7CB and similar dimers over a number of years, going back to the original synthetic work where an additional phase of unknown structure was identified below a conventional nematic. It was later suggested that this might correspond to the ${\mathrm{N}}_{\mathrm{TB}}$ phase predicted by Meyer [38] and by Dozov [39]. Considerable experimental and modelling work (including atomistic simulations [40]) has now taken place, confirming, almost without doubt, that CB7CB shows a chiral phase structure with domains of opposite handedness [40,41,42,43] despite the fact that the molecules themselves are achiral. Up to a few years ago, atomistic simulation models were not possible for this type of phase, but good simulation models are now able to provide a molecular level picture of the order within such phases, subject to the usual constraints of system size and simulation time scale. Simulations confirm the transition from a nematic to a chiral twisted structure associated with an order parameter change, as shown in Figure 4, and confirm a helical pitch, order parameter, and conical tilt angle in agreement with experiment [44].

In recent work, Boyd and co-workers used atomistic simulations to study the local ordering of some very unusual liquid crystal phases where the structure of the phase ultimately arises from local molecular packing. These studies include the following:

• The dark conglomerate (DC) phase where atomistic simulations demonstrate the presence of molecular layers which undergo a saddle-splay layer deformation [31];
• Studies of B4 phase forming molecules [45,46] where subtle changes in molecular packing can be induced by the positioning of a (sterically significant) lateral methyl group at chiral centres in bent-core molecules with a long and short arm. These lead to changes in the splay, twist, and bend of molecular layers, which ultimately lead to the formation of twisted filaments and multi-level hierarchical self-assembled structures.

In the early days of liquid crystal simulation, simple coarse-grained models such as hard and soft spherocylinders, hard ellipsoids, the Gay–Berne model [11,12,13,47,48,49], and models composed of joined hard spheres [50,51] were instrumental in providing insights into simple nematic and smectic phases and polymer liquid crystals [20,22,52,53]. It is an interesting question whether a new generation of models can provide further insights into more complex phases, where phase structure appears to be dominated by local packing effects. Already considerable success has been achieved by exploring the influence of packing effects in non-linear rigid particles in a number of key areas:

• The formation of chiral superstructures from asymmetric bent-cores molecules composed of achiral tangential Lennard–Jones and WCA spheres [54];
• The formation of the ${\mathrm{N}}_{\mathrm{TB}}$ phase from crescent-shaped particles composed of tangential hard spheres [55];
• The formation of biaxial, twist-bend, and splay-bend nematic phases from hard banana-shaped particles [56].

The latter work is particularly interesting because it shows that in this case biaxial, twist-bend, and splay-bend nematic phases are metastable with respect to smectic phases but stability can be induced by the presence of polydispersity in the particle’s length or by curvature in the particle shape. This is immediately relevant to the design of new colloidal liquid crystals [57] but is also pertinent to understanding many thermotropic systems; i.e., many of the new phases that exhibit chiral arrangements of molecules form in systems where shape dominates the packing of molecules but where there is some “polydispersity” in molecular shape arising from conformational disorder that destabilises both crystal and smectic phase structures.

Figure 5 shows preliminary work from the work of the Wilson group exploring new coarse-grained models. Figure 5a shows a snapshot from a model composed of three spherocylinders interacting through the soft spherocylinder model of Lintuvuori and Wilson [23,58]. Here, some “polydispersity” in shape arises from flexibility built into the model through a torsional potential about the central spherocylinder. Figure 5b shows a snapshot from a rigid model of a chiral asymmetric-bent-core molecule. In this model, microphase separation is ensured through the use of Lennard–Jones sites to represent the central aromatic part of the molecule and WCA sites to represent aliphatic parts of the molecule, leading to smectic layer formation. Packing frustration is ensured by requiring a larger volume to pack aliphatic tails and through the presence of a chiral off-axis site that mimics the effects of lateral methyl substitution in the work of Hegmann and co-workers [45,46]. Here, twisted and splayed smectic layers are induced from these packing constraints, leading to a helical superstructure formed from twisted layers forming throughout the phase.

2.3. Calculation of Material Properties Using Molecular Simulations

In addition to research aimed at predicting the structure of liquid crystal phases and their transition temperatures, considerable effort has been invested into developing methods that might be suitable for predicting material properties. The latter is potentially important for future screening applications, where simulation could act as a valuable tool in helping design appropriate liquid crystal molecules prior to synthesis based on predictions of their material properties. Much of the initial work in this area has been carried out on coarse-grained models, which can be simulated using large numbers of molecules and studied over long simulation times. Both of these criteria tend to be important for reliable calculations of a large number of key material properties.

Of particular note has been the work carried out in developing predictive methods for Frank elastic constants [59,60,61]. Here, the most successful approach [59] considers wavevector-dependent fluctuations in orientational order, which occur naturally during the course of a molecular dynamics simulation of a nematic liquid crystal at finite temperature. This approach has been successfully applied to nematic phases composed of hard prolate ellipsoids and hard spherocylinders [61,62] and to Gay–Berne mesogens [59,60]. However, the large system sizes required for these calculations make this methodology very challenging to implement for atomistic models at the current time. In principle, elastic constants can also be obtained from an approximate approach involving the direct simulation of the Freédericksz transition [63,64] (although this has proved to be inefficient compared to the fluctuation approach for simple systems and would be very challenging to carry out accurately for atomistic studies) or from calculations of the direct correlation function in the nematic phase [65].

Some progress has also been made in the calculation of other display-related material properties. An example is the rotational viscosity coefficient, ${\gamma }_1$, which is important (for example) in determining the on and off times of a twisted nematic display ($t_{\mathrm{off}}\propto {\gamma }_1$, $t_{\mathrm{on}}\propto {\gamma }_1$) and for other nematic devices. ${\gamma }_1$ can be determined by a number of techniques including equilibrium and non-equilibrium molecular dynamics methods [66,67,68,69]. While extensive work has been carried out on simple coarse-grained models, such as Gay–Berne potentials, these simulated values do not closely match the values of ${\gamma }_1$ measured for real thermotropic systems [69]. To some extent, this is because most Gay-Berne parametrisations more closely match colloidal liquid crystals than small thermotropic organic molecules. However, equilibrium molecular dynamics can be applied successfully, and relatively easily, to atomistic studies [16], and have the advantage of avoiding the need to couple molecules to an external field. Here, for example, ${\gamma }_1$ can be successfully obtained from the director angular velocity correlation or the director mean-squared displacement. Both methods rely on thermally excited fluctuations in orientational order, avoiding the need to perturb the system by application of an external field.

Atomistic studies have also been successful in calculating flexoelectric coefficients for a nematic liquid crystal [17]. Here, equilibrium molecular dynamics calculations have been performed for a series of temperatures in the nematic phase. The simulation data are used to calculate the flexoelectric coefficients $e_{\mathrm{s}}$ and $e_b$ using the linear response formalism of Osipov and Nemtsov [70].

A further interesting area where atomistic simulations have been applied successfully is in the study of helical twisting powers. The helical twisting power (HTP) measures the ability of a chiral molecule to impart a chiral twist to a bulk nematic phase. Two successful approaches have been demonstrated involving the chirality order parameter approach of Ferrarini and co-workers [71,72,73,74,75,76,77] and a scaled chirality index arising from the work of Osipov and co-workers [78] that has been employed for molecular systems [79,80]. Part of the success of the two methods arises from them being single-molecule techniques; i.e., they do not require the simulation of bulk phases. Both approaches have been extensively studied and the results verified against experiments for a range of chiral molecules [81,82].

Interestingly, the chirality order parameter approach has also been applied to achiral bent-core molecules [83]. It is noted that some achiral bent-core molecules, acting as dopants, can increase the helical twisting power of a chiral host phase, i.e., the exact opposite of what would be expected. Calculations of the chirality order parameter, $\chi$, show that although these molecules are on average achiral, they can exhibit “chiral conformations, which have extremely large values of $\chi$ in comparison to standard chiral dopants [83,84,85]. Hence, the suggestion is that these conformations have such large twists that they lead to a spontaneous asymmetry between right and left-handed forms in a chiral environment, through chiral templating with solvent molecules or other bent-core molecules (i.e., the handedness of the phase leads to there being a slight preference for either left-handed or right-handed conformations for the bent-core molecule). This then increases the overall chirality of the bulk phase [83]. This hypothesis was tested in bulk simulations of a simple coarse-grained model that exhibited high helical twisting powers for chiral conformations [86].

In principle, the HTP can also be deduced from the sampling of torques exerted on neighbouring molecules in a bulk phase [87,88]. This is a far harder calculation than a single-molecule based approach to calculating HTP, as it involves the simulation of bulk phases and the calculation of torques. Its application to atomistic models has been limited to rigid molecules within a coarse-grained Gay-Berne nematic solvent but the method does yield the correct sign of the helical twist induced by a molecule and provides a good approximation to the magnitude of HTP [87].

It is also possible to look at individual chiral dopant molecules within a uniformly twisted nematic phase of wavevector $k=2\pi /P$ (where P is the pitch of the twisted phase). Here, there should be a chemical potential difference between enantiomers that want to induce a twist in the same direction as the host phase and those that want to induce a twist in the opposing direction. In the limit of infinite dilution this difference in chemical potential, $\Delta \mu$, is directly proportional to the helical twisting power, $\beta$, and the twist elastic constant $K_2$ [89].

$$\Delta \mu ={\mu }_{-}-{\mu }_{+}=8\pi \beta K_{2}k$$

(4)

This was exploited by Wilson and Earl using twisted periodic boundary conditions to calculate the helical twist power of five atomistic dopant molecules in simple coarse-grained solvents [90]. Relatively large errors arise with this method, as the dopants must be “grown” into the solvent in a series of separate simulation steps. In the original paper, this was performed for spherocylinder and Gay–Berne solvents. However, for a uniformly twisted nematic solvent, the key quantity of interest is the twist elastic constant, and the exact details of the solvent potential are not important, providing (of course) a uniform nematic is simulated. Therefore, in principle, soft-core model potentials [58] could be used to model the liquid crystal solvent, greatly reducing the difficulty of introducing a chiral dopant into a bulk phase in which there could initially be strong overlaps between the dopant and the solvent molecules.

Finally, it is worth noting that only very recently have advances in computer power and parallel algorithms opened up the possibility of studying large atomistic simulations of several thousand molecules. As computer time continues to increase in the future (and atomistic simulations increase further in size), it should be possible to look at other material properties that are currently difficult to study because of system size limitations. Dielectric anisotropy and ion conductivity in liquid crystals are good examples of such properties.

3. Lyotropic Liquid Crystals

3.1. Surfactant Models, Micelles, and the Formation of Lyotropic Phases

Conventional lyotropic liquid crystals have traditionally been very difficult to study using molecular simulation methods. At an atomistic level, the self-assembly of surfactants to give micelles and, subsequently, the self-assembly of micellar aggregates to give liquid crystal phases, occurs on time scales that are extremely challenging for atomistic simulation. However, advances in computer time mean that it is now possible to “see” the self-assembly process occurring for conventional amphiphiles for small aggregates if calculations are carried out well above the critical micelle concentration (CMC) and simulations are extended for tens, or (in some cases) hundreds, of nanoseconds. Moreover, coarse-grained simulation methods have become very useful in studying traditional lyotropic phase diagrams. Figure 6 shows a micelle that has spontaneously self-assembled from a dispersed group of monomers for the industrially important surfactant LAS (linear alkylbenzene sulfonate). Here, simulations were carried out at the fully atomistic level using the GAFF force field in combination with TIP3P water using the completely linear isomer of LAS shown in Figure 6a (noting that in most commercial applications LAS is used as a mixture of branched isomers a shown in Figure 6b). Here, micelles form within 100 ns of simulation for a system with a LAS:water ratio of 1:221, which is well above the CMC.

For many surfactant systems, the CMC occurs at sufficiently low concentrations that atomistic studies are not feasible near the CMC because of the number of water molecules required and also the time scales (beyond the microsecond regime) required for diffusion of water molecules. However, good progress has been made using coarse-grained models. Figure 6 additionally shows models and snapshots from a typical dissipative particle dynamics (DPD) coarse-grained surfactant study. In DPD, soft spheres are used to represent coarse-grained sites corresponding to groups of atoms within the molecule. Here, soft spheres interact through a conservative force:

$$\begin{array}{ccc}\hfill F_{\mathrm{C},ij}& =& a_{ij}(1-{\mathbf{r}}_{ij}/r_{\mathrm{cut}})\mathrm{for}\left|{\mathbf{r}}_{ij}\right|<r_{\mathrm{cut}}\hfill \\ \hfill F_{\mathrm{C},ij}& =& 0\mathrm{for}|{\mathbf{r}}_{ij}|\ge r_{\mathrm{cut}}\hfill \end{array}$$

(5)

which leads to a harmonic repulsive potential. The chemical interactions are mainly controlled by the parameters $a_{ij}$, which govern the strength of the repulsive interactions between two species i and j. Part of the beauty of DPD is that simulations can be carried out using large dynamic time steps (due to the use of soft spheres and soft springs linking them) with a greatly reduced number of sites in comparison to atomistic simulation. This comes with some loss in terms of chemical specificity. However, in recent years DPD has become a powerful methodology for modelling surfactant systems and considerable work has been carried out in mapping the interactions in DPD models to “real” molecular systems [92,93,94,95,96].

Figure 6d shows three models for different LAS molecules that provide good phase diagram predictions for linear and branched LAS species [91]. Figure 6e shows snapshots taken from DPD simulations of a linear LAS model in the high-concentration regime, which corresponds to hexagonal phases composed of rod-shaped micelles and (at slightly higher concentrations) lamellar phases. The lamellar layers exhibit a spacing of between 3.1 and 3.3 nm, in excellent agreement with experiments (3.2–3.3 nm).

3.2. Models for Chromonic Liquid Crystals

Of significant interest over the last decade has been the study of chromonic lyotropic liquid crystals [97]. These systems arise from non-conventional amphiphiles in which solubilizing groups are attached to the periphery of “molecular discs” (see Figure 7). Chromonic mesogens have been seen in a variety of dye and drug molecules, which commonly exhibit these two structural motifs (discs and hydrophilic peripheral groups). Chromonics have recently found new potentials applications [98] in areas such as biosensors [99,100], in the fabrication of thin-film structures [101] and in controllable self-assembly of gold nanorods [102], and interest in chromonic has been extended to drug delivery systems and to new areas such as “living liquid crystals” [103]. Understanding chromonic self-assembly is important in all these areas.

For chromonics systems, simple coarse-grained models have demonstrated the formation of the most commonly seen mesophases: the chromonic N and M phases. Simulations have confirmed the structure of these mesophases and have allowed for the mapping out of chromonic phase diagrams as a function of concentration [104,105].

Alongside coarse-grained studies, atomistic simulations have been productive in the following:

• Clarifying the nature of self-assembled aggregates in dilute solution [106];
• Explaining the origin of different types of stacking motifs [107];
• Explaining the formation of “layer aggregates”, which can lead to smectic or layered chromonic phases [107,108].

These are phenomena that are difficult to probe experimentally at a molecular level in dilute solutions and for which simulation thereby provides unique information.

Figure 7. (a) The disc-shaped structure of the chromonic dianionic monoazo food dye sunset yellow in the NH hydrazone tautomer (the stable tautomeric form in aqueous solution); (b) stacking of sunset yellow molecules in solution into a chromonic aggregate (blue spheres represent sodium counterions) as seen in molecular dynamics simulations of sunset yellow [106]; (c) schematic diagram showing typical alignment of chromonic aggregates to form a liquid crystalline phase.

Figure 7. (a) The disc-shaped structure of the chromonic dianionic monoazo food dye sunset yellow in the NH hydrazone tautomer (the stable tautomeric form in aqueous solution); (b) stacking of sunset yellow molecules in solution into a chromonic aggregate (blue spheres represent sodium counterions) as seen in molecular dynamics simulations of sunset yellow [106]; (c) schematic diagram showing typical alignment of chromonic aggregates to form a liquid crystalline phase.

3.3. Studies of Nonionic Chromonics

The structure shown in Figure 8a is of 2,3,6,7,10,11-hexa-(1,4,7-trioxa-octyl)-triphenylene (TP6EO2M), which consists of a central polyaromatic core (a triphenylene ring) functionalised by six hydrophilic ethyleneoxy (EO) chains. TP6EO2M is the archetypal nonionic chromonic where the hydrophobic interactions arising from the aromatic rings induce stacking into a chromonic column but the hydrophilic interactions of the short EO chains are sufficient to solubilise the resulting aggregates.

TP6EO2M has become the “fruit fly” of chromonic liquid crystals simulations and has been simulated by a range of models, as shown in Figure 8. This is partly because the structure is relatively simple and symmetrical but also because the formation of chromonic stacks and chromonic phases for this molecule provides a major challenge for modern methods of multi-scale modelling.

All-atom simulations of TP6EO2M (Figure 8a) using the OPLS-AA force field [109] demonstrate the formation of chromonic stacks. The free energy for association, $\Delta G_{\mathrm{agg}}^{°}$, can be approximated from the depth of the “attractive well” obtained from potential of mean force (PMF) calculations, $U_{\mathrm{PMF}}\left(r\right)$. Here, a molecule is pulled away from a dimer, or more generally from a n-mer, to a point where they are no longer interacting. In practice these can be performed by numerically integrating the average constraint force, $f_{\mathrm{c}}$, obtained at a series of separation distances, s:

$$U_{\mathrm{PMF}}\left(r\right)={\int }_r^{r_{\max }}\left[{〈f_{\mathrm{c}}〉}_s+\frac{2k_{\mathrm{B}}T}{s}\right]ds$$

(6)

where r is the distance for the PMF and $2k_{\mathrm{B}}T/s$ is a kinetic entropy term, which accounts for the increase in rotational volume at larger separation distances [110,111,112]. Many chromonic mesogens are assumed to undergo isodesmic association, with the same free energy change for each subsequent addition of a molecule to an n-mer. This leads to an exponential distribution of aggregate sizes. For TP6EO2M, the binding energy for a n-mer ($\Delta G_{agg}^{°}=-12RT$) is in good agreement with experiments. However, for TP6EO2M, and for many other chromonics, we find “quasi-isodesmic” association where the binding free energy, $\Delta G_{\mathrm{agg}}^{°}$, for two molecules to form a dimer is slightly larger in magnitude (∼−2.5 $RT$ lower for TP6EO2M) than the binding energy for an n-mer. Typically, for molecules in the interior of a stack, the entropy loss from chain confinement and orientational ordering upon aggregation is larger than on the formation of a dimer when the molecules are only confined by the presence of one neighbour.

Values of the enthalpy and the entropy of association can both be obtained from the temperature dependence of $\Delta G_{\mathrm{agg}}^{°}$.

$$\Delta H^{°}=R\left[\frac{\partial \left(\Delta G_{\mathrm{agg}}^{°}/RT\right)}{\partial \left(1/T\right)}\right].$$

(7)

Interestingly, for TP6EO2M, $-T\Delta S^{°}$ and $\Delta H^{°}$ both contribute favourably to $\Delta G_{\mathrm{agg}}^{°}$ and at 300 K it is the entropic contribution that is larger, indicating that the confinement of chains and orientations on molecular association is more than compensated for by the gain in entropy of the solvent that is released from interacting with the phenyl rings of triphenylene. This finding is very much in-line with traditional interpretations of the hydrophobic effect for small molecules at room temperature.

As might be expected, the relative balance of hydrophilic and hydrophobic interactions in these atomistic simulations is critical, and “out-of-the-box”, the GAFF force field which normally performs well for a range of ionic chromonic systems behaves badly in comparison to OPLS, with the formation of disordered globular aggregates. With GAFF, the EO units are insufficiently hydrated to stabilise chromonic stacks.

At the other end of the scale of models in Figure 8c is a dissipative particle dynamics (DPD) representation of TP6EO2M: A model was designed by Walker et al. [104]. For TP6EO2M, two things are critical: the balance of interactions between hydrophobic and hydrophilic beads and also $\Delta G_{\mathrm{agg}}^{°}$ relative to the DPD energy scale, i.e., $ϵ$ that controls the magnitude of the reduced temperature $k_{\mathrm{B}}T/ϵ$. The former can be tuned for many DPD systems based on experimental infinite dilution activities or chemical potentials obtained from dilute atomistic simulations. For TP6EO2M, the latter is available from atomistic modelling (as discussed above).

DPD is sufficiently tractable to allow the full chromonic phase diagram to be simulated for TP6EO2M, producing an approximately exponential distribution of stacks in the isotropic phase (in-line with quasi-isodesmic association); a nematic N phase at intermediate concentrations, where columns are sufficiently long to align due to excluded volume interactions favouring translation entropy over orientational entropy; and a hexagonal M phase at higher concentrations, where the stacks pack on a hexagonally ordered lattice.

If EO chains are replaced in this model with hydrophobic–lipophobic chains [105] (e.g., as might occur, for example, with siloxanes or fluorinated systems) complex supramolecular aggregates form in which hydrophobic–lipophobic chains are excluded from water by the joining together of stacks with a single-molecule cross-section into dimers or trimers. This provides a mechanism for the formation of a novel chiral aggregate and also in the case of “Janus mesogens” with three adjacent hydrophobic–lipophobic chains, a novel smectic chromonic phase.

It is also interesting to ask whether an “in-between” coarse-grained model can provide further insights into the behaviour of TP6EO2M. This was the aim of recent work by Potter, Wilson, and co-workers who explored a MARTINI-style coarse-grained model (Figure 8b) where sites are grouped into beads representing two or three heavy atoms plus attached hydrogen atoms [113,114]. Such models are quite interesting as they aim to capture the individual interactions responsible for local order but also aim to capture the thermodynamics and correct self-assembly and self-organisation over larger length scales. They aim to perform this while also hitting a “sweet spot” where they employ only a fraction of the simulation time of an atomistic model, saving CPU cycles through a reduced number of sites, long dynamic time steps, and a speed-up in phase space exploration. While the idea of such a “Goldilocks” model is very appealing and has also many potential uses in related fields such as polymer simulation, membrane studies [5,7,115,116,117,118,119] and protein simulation, the challenges of making such a model are quite demanding. In fact, the very act of coarse-graining in an aqueous system alters the balance between entropy and enthalpy. The hydrophobic interaction in coarse-grained models of chromonics must incorporate additional enthalpic contributions to the free energy to make up for the loss of entropic ones.

There are a range of methods available to generate an intermediate coarse-grained model for a molecule such as TP6EO2M using both “bottom-up” and “top-down” routes [120], i.e., starting from an atomistic reference simulation, or starting from experimental structural and thermodynamic data (see the diagram in Figure 9).

Traditionally, structure-based coarse-graining by methods such as multiscale coarse-graining (MS-CG) in the form of hybrid force matching (HFM) [120,121,122,123] or using iterative Boltzmann inversion [124], do a good job of capturing “molecular structure” and “order” but produce coarse-grained potentials that do not capture thermodynamics well (particularly mixing thermodynamics). They also tend to produce CG potentials that are extremely density and temperature-dependent and so have very limited transferability between state points. While HFM for TP6EO2M using an atomistic reference stack [113] yields quite good structures for the CG model, the very high binding energies make this rather unsuitable as a method to study changes in the aggregation of chromonics molecules with concentration.

It is worth noting that successful chemical engineering theories for mixtures, such as SAFT (Statistical Associating Fluid Theory), are often useful in predicting phase diagrams for simple molecular mixtures. These are starting to become useful in providing “top-down” pathways to thermodynamically consistent coarse-grained models. Of particular interest here is the SAFT-$\gamma$ Mie model [2,125,126,127,128,129,130], which represents coarse-grained sites by generalised Mie potentials:

$$U_{\mathrm{Mie},ij}=C{ϵ}_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r}\right)}^{{\lambda }_{\mathrm{r}}}-{\left(\frac{{\sigma }_{ij}}{r}\right)}^{{\lambda }_{\mathrm{a}}}\right]$$

(8)

where exponents ${\lambda }_{\mathrm{r}}$ and ${\lambda }_{\mathrm{a}}$ govern the ranges of the repulsion and attraction of the potential, and the following is the case.

$${ϵ}_{ij}=\left(1-k_{ij}\right)\frac{\sqrt{{\sigma }_{ii}^3{\sigma }_{jj}^3}}{{\sigma }_{ij}^3}\sqrt{{ϵ}_{ii}{ϵ}_{jj}}.$$

(9)

For species that are chemically incompatible for entropic or enthalpic reasons, $k_{ij}>0$, and for species where one group is well solvated by another, $k_{ij}<0$ (recalling that within a coarse-grained model an entropically unfavourable interaction is often partially translated into an unfavourable enthalpic interaction due to the reduction in the degrees of freedom). SAFT was exploited by Potter et al. to develop a SAFT-$\gamma$ Mie model for TP6EO2M [114]. The balance of $k_{ij}$ values was found to be crucial. If ethylene oxide units are poorly solvated, chromonic stacks do not form and instead phase separation, or a conglomerate of aggregates, occurs. Moreover, $k_{ij}$ must be greater than zero for both aromatic-water interactions and for aromatic-ethylene oxide interactions in order to see stable chromonic stacks. In other cases, no aggregation occurs, or the system phase separates with the formation of large non-structured aggregates depending on specific $k_{ij}$ combinations. Chromonic aggregation, therefore, appears only in a small region of parameter space where the hydrophilic–hydrophobic balance between aromatic, ethylene oxide, and water is just right.

It is perhaps not surprising that the most successful intermediate coarse-grained model for TP6EO2M comes from a Martini 3 model of the mesogen [113]. Martini models have been parametrised to reproduce good solvation free energies for different coarse-grained sites and, hence, are able to capture the correct hydrophilic–hydrophobic balance required to see chromonic behaviour. Figure 10 shows two snapshots from simulations of TP6EO2M using the tuned Martini 3 model of Potter [113] and co-workers, showing the structure of chromonic N and M phases. The latter occurs at high mesogen concentrations when, as shown in Figure 10, the columns are sufficiently close to undergo hexagonal packing.

3.4. Ionic Chromonics: A Rich Variety of Aggregation Motifs

Yu and Wilson [107] have used atomistic models to explore the rich variety of chromonic aggregation that occurs in cyanine dye systems, studying four dyes: pseudoisocyanine chloride (PIC), pinacyanol chloride (PCYN), 5,5${}^{\prime }$,6,6${}^{\prime }$-tetrachloro-1,1${}^{\prime }$,3,3${}^{\prime }$-tetraethylbenzimidazolyl-carbocyanine chloride (TTBC), and 1,1${}^{\prime }$-disulfopropyl-3,3${}^{\prime }$-diethyl-5,5${}^{\prime }$,6,6${}^{\prime }$-tetrachloro-benzimidazolylcarbocyanine sodium salt (BIC). Simulations with the GAFF/TIP3P force field combination allowed a range of aggregate behaviour to be observed and rationalised based on potentials of mean force obtained from pulling molecules from the end of stacks of two, three, or four molecules. Both H-aggregates (where molecules exhibit face-to-face stacking and show hypsochromic (blue) spectra shifts in comparison to monomers) and J-aggregates (where adjacent molecules exhibit staggered stacking and show bathochromic (red) spectra shifts) are seen depending on the system, with molecular stacking arrangements leading to Y-junction and shift defects and sheet-like assemblies of molecules.

Examples of the aggregate structures seen are shown in Figure 11. What is particularly interesting is the layered structures seen by Yu and Wilson and also by Thind and co-workers [108], which correspond to a different type of chromonic aggregate to that normally reported. This layered structure can form the basis for a chromonic smectic phase. Moreover, it is likely that the layered structures reported represent an intermediate state that precedes the organisation of molecules into tubular or other higher-order architectures. In experimental studies, large scale tubular structures have been reported for some cyanine dyes using electron microscopy techniques [131,132,133].

It is extremely challenging to develop coarse-grained molecular models that are able to capture the same types of aggregation seen with PIC, PCYN, TTBC, and BIC, including the thermodynamic equilibrium between monomers, H-aggregates and J-aggregates: M ⇌ H ⇌ J. However, such models are required to explore larger-scale aggregation beyond the time and length scales that are accessible to atomistic studies.

Yu and Wilson have explored this for one perylene dye molecule bis-(N,N-diethylaminoethyl)perylene-3,4,9,10-tetracarboxylic diimide dihydrochloride, PER (Figure 12) [8]. While bottom-up models fail to capture the correct behaviour of PER molecules in solution, it is possible to tune a top-down Martini 3 model to capture the same aggregation behaviour seen by atomistic studies. Moreover, a proof-of-concept hybrid model combining hybrid force matching (HFM) and Martini 3 was successfully developed. Here, a combined approach is able to capture the aggregate structure and provide a good representation of short-range electrostatic interactions that often control the local structure of chromonic aggregates. The use of Martini water and ions provides a large reduction in the number of sites within the model and improves the thermodynamics of hydration (which is usually poorly represented within a force matching approach). In Figure 12, we show the atomistic and coarse-grained structures seen by Yu and Wilson, together with the potential of mean force corresponding to the binding of dimer molecules. The phase structures shown in the lower part of Figure 12 are seeded to produce columns that connect over the periodic boundary conditions, i.e., effectively columns of infinite length. This makes it possible to investigate the concentration corresponding to the N to M phase transition, when the concentration is sufficiently high to see hexagonal packing of columns.

4. Conclusions

The huge increase in computer time that has become available over the last decade, together with new algorithms capable of exploiting parallel computer architectures, has opened up many possibilities for the modelling of both thermotropic and lyotropic liquid crystal systems. This perspective article, which is part of a collection of articles showcasing developments within the UK, has pointed to several exciting areas of research emerging.

For thermotropic nematics, it is now possible to simulate a few thousand molecules, and, with the advent of improved molecular force fields, provide predictions for nematic clearing temperatures, ${\mathrm{T}}_{\mathrm{NI}}$ within a 10 ${}^{°}$C (or better) range (see Section 2.1). Current studies are limited to a few hundred nanoseconds of simulation time for this size of system. Nonetheless, as shown in Section 2.2, this type of simulation is also able to provide valuable insights into the structure of more complex liquid crystalline phases, where the details of molecular arrangements may be very difficult to probe experimentally. Section 2.2 also points to the possibilities of using simpler coarse-grained models to study larger systems of molecules as a method of understanding the structure of complex layered phases and/or thermotropic phases where chirality arises from achiral systems. Here, packing arrangements between molecules are often very important but there are many challenges arising from the complex interplay of packing constraints, microphase separation and favourable local interactions arising from the presence of specific functional groups.

Considerable progress has now been made in the development of methods for the calculation of liquid crystalline material properties for thermotropics. Undoubtedly, coarse-grained models have proved extremely useful as a testbed for this work, and it has been possible to test the convergence of results over long simulations with large simulation sizes. It should be noted that some of these studies point to the need to carry out simulations that are larger than can reasonably be simulated at an atomistic resolution. However, as computer time increases, this restriction will gradually be removed. Already good progress has been made in the calculation of helical twisting powers, rotational viscosities, and flexoelectric coefficients at an atomistic level.

For conventional lyotropic systems, formed through surfactants, atomistic simulations are often limited to the study of single micelles, which can be seen to form spontaneously at concentrations well above the CMC. However, recent years have seen many developments in coarse-grained models for these systems, particularly in the area of dissipative particle dynamics. These models are able to provide a good link to underlying thermodynamic data measured experimentally or determined by theory or modelling. They can also provide good predictions for experimentally measured observables, such as the lamellar layer spacing (see Section 3.1).

Chromonic liquid crystals were also discussed in detail. For these, non-conventional amphiphiles atomistic simulations have been particularly productive in clarifying the nature of the self-assembled aggregates in dilute solution (Section 3.2), which have the potential to be more complex than those seen for conventional surfactants. In recent years, different forms of self-assembly have been noted particularly in chromonic dye molecules, from the most common form of association involving stacking into columns (H-aggregate formation), to defect stacking with shift-defects and Y-defects (leading to branched structures), and self-assembly into brickwork layers (Section 3.4).

Considerable progress has also been made in the development of different levels of coarse-grained models for chromonic liquid crystals (see Section 3.3 and Section 3.4). Of significant interest here is the challenge of capturing the correct local stacking arrangements at the same time as capturing the correct thermodynamics of self-assembly. Noting that in coarse-graining chromonic systems, the balance of entropy and enthalpy shifts as degrees of freedom are removed from the system. Coarse-grained models at DPD and bead-spring level have both proved successful in visualising the chromonic N and M phases.

This last decade has seen the modelling of liquid crystalline systems mature considerably. The decade ahead promises much in terms of developments in modelling. For liquid crystals, machine learning methods are likely to further improve molecular force fields. Moreover, such methods are likely to lead to a new generation of coarse-grained models. Better connection of models across time and lengths scales may well lead to truly multi-scale modelling approaches for understanding liquid crystal systems.