crystals

Article
Molecular Simulation Approaches to the Study of Thermotropic
and Lyotropic Liquid Crystals
Mark R. Wilson * , Gary Yu , Thomas D. Potter, Martin Walker, Sarah J. Gray, Jing Li and Nicola Jane Boyd

Chemistry Department, Durham University, Durham DH1 3LE, UK; [Link]@[Link] (G.Y.);
[Link]@[Link] (T.D.P.); [Link]@[Link] (M.W.); [Link]@[Link] (S.J.G.);
[Link]@[Link] (J.L.); [Link]@[Link] (N.J.B.)
* Correspondence: [Link]@[Link]

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.
Citation: Wilson, M.R.; Yu, G.; Potter,
The development of effective coarse-grained models for chromonics is discussed in detail, including
T.D.; Walker, M.; Gray, S.J.; Li, J.;
models that have demonstrated the formation of the chromonic N and M phases.
Boyd, N.J. Molecular Simulation
Approaches to the Study of
Thermotropic and Lyotropic Liquid
Keywords: liquid crystals; molecular simulation; molecular dynamics; dissipative particle dynamics
Crystals. Crystals 2022, 12, 685.
[Link]
cryst12050685
1. Introduction
Academic Editor: Ingo Dierking
In recent years, molecular simulation has become a powerful tool for studying a range
Received: 22 March 2022 of soft matter systems. Simulations have been extended to bulk polymers [1], polymer
Accepted: 30 April 2022 surfaces [2], proteins [3,4], membranes [5–7], self-assembly in solution [8], and molecules
Published: 10 May 2022 at water interfaces [9,10], in addition to many other systems. In these cases, simulations
Publisher’s Note: MDPI stays neutral aim to provide quantitative predictions to compare with experiment and also to provide
with regard to jurisdictional claims in qualitative insights into local molecular structure and order, which are often difficult to
published maps and institutional affil- obtain by experimental means. These comments are particularly true for liquid crystalline
iations. 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–14].
Molecular simulation models are part of a traditional hierarchy of simulation method-
Copyright: © 2022 by the authors. ologies which cover different time and length scales from the microscopic to mesoscopic to
Licensee MDPI, Basel, Switzerland. continuum scales. For liquid crystals, this hierarchy is particularly significant (see Figure 1)
This article is an open access article
because it covers the following:
distributed under the terms and
conditions of the Creative Commons
• The quantum mechanical regime, where single-molecule calculations are valuable in
Attribution (CC BY) license (https://
determining molecular properties of single thermotropic mesogens to use, for example,
[Link]/licenses/by/ in the development of materials for displays [15];
4.0/).

Crystals 2022, 12, 685. [Link] [Link]

Crystals 2022, 12, 685 2 of 24

• 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–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].
macroscopic

Continuum
Models
Lattice
mesoscopic
length

models
Coarse-grained
simulation
microscopic

Atomistic
simulation
Quantum
mechanical

ps ns μs ms s
time
Figure 1. The time and length scale hierarchy for liquid crystal simulations covering microscopic,
mesoscopic, and macroscopic regimes.

Over the last decade, the possibilities for modelling liquid crystals have been extended
hugely with the realisation of readily accessible parallel computer architectures and sim-
ulation 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 develop-
ments 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 tem-
peratures, 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!
Crystals 2022, 12, 685 3 of 24

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–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:

EMM = ∑ Kr (r − req )2 + ∑ Kθ (θ − θeq )2

bonds angles
5
+ ∑ ∑ Cn (cos(ψ))n + ∑ k d (1 + cos(nd ω − ωd )) (1)
torsions n=0 impropers
N   12 6  
1 qi q j
σ
∑
σij ij
+ 4ε ij rij − rij + 4πε0 rij ,
i>j

where req and θeq are, respectively, natural bond lengths and angles, Kr , Kθ , and Cn are,
respectively, bond, angle, and torsional force constants, ψ is a dihedral (torsional) angle, k d
is a force constant associated with improper dihedral angles, ω, nd , and ωd , respectively,
represent an improper dihedral angle, its periodicity, and the phase angle associated with
it. σij and ε ij are the usual Lennard–Jones parameters modelling nonbonded interactions,
and qi , q j are partial electronic charges. Here, we have written EMM 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 dramati-
cally 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–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
Crystals 2022, 12, 685 4 of 24

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.

Figure 2. Simulation snapshots of 2048 molecules from the nematic phase of the bent-core mesogen
C5-Ph-ODBP-Ph-OC12 at 480 K. (Left): Line drawing representation of the molecular bent core within
the nematic phase. (Right): space-filling representation of C5-Ph-ODBP-Ph-OC12 molecules in the
nematic phase showing molecular cores in green and alkyl tails in gold. The snapshot shows the
beginnings of microphase separation between cores and tails that occurs in a pretransitional region
before the phase transition to a DC phase at lower temperatures.

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 ~n and an orientational order
parameter S2 . In practice, these are obtained by calculating the ordering tensor:

N
1
∑
 
Qαβ (t) = 3uiα uiβ − δαβ , α, β = x, y, z, (2)
2N i =1

where the sum runs over all N molecules. The largest eigenvalue of the Q tensor represents
the following:
1 N
N i∑
P2 (t) = P2 (cos θi ), (3)
=1
where P2 is the second Legendre polynomial, and the associated eigenvector is the director
~n(t). S2 is defined as the time average of P2 over a suitable time interval where h P2 i is
unchanging. However, in practice, to minimise system size effects in locating the phase
transitions, S2 is often obtained from −2× and the middle eigenvalue of Q, which fluctuates
about a value of zero in the isotropic phase but equals P2 (t) 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.
Crystals 2022, 12, 685 5 of 24

For the systems tested (molecular structures given in Figure 3), GAFF-LCFF provides
excellent liquid crystal clearing point predictions, as shown in Table 1. 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].

Figure 3. Structures of the five mesogens used in Table 1 (top to bottom: 1,3-benzenedicarboxylic
acid,1,3-bis(4-butylphenyl)-ester, C5-Ph-ODBP-Ph-OC12, C4-Ph-ODBP-Ph-C7, C4O-Ph-ODBP, and
C4O-Ph-ODBP (trimethylated)).

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].

TNI (exp.) TNI (GAFF-LCFF)

Molecule
/K /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

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 experi-
ments. Figure 4 is from the preliminary work of Yu and Wilson and shows the structure
of the twist-bend nematic (NTB ) 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
Crystals 2022, 12, 685 6 of 24

below a conventional nematic. It was later suggested that this might correspond to the
NTB 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–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].

Figure 4. (Left): Snapshot, with orientational colour coding, showing the structure of the NTB phase
of CB7CB. (Right): The structure of a CB7CB molecule.

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–13,47–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:
Crystals 2022, 12, 685 7 of 24

• The formation of chiral superstructures from asymmetric bent-cores molecules com-

posed of achiral tangential Lennard–Jones and WCA spheres [54];
• The formation of the NTB 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 sphero-
cylinders 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.

(a) (b)
Figure 5. (a): Snapshot from a simulation of a simple coarse-grained simulation model of jointed
spherocylinders. 1728 molecules are used with each molecule composed of three bonded sphero-
cylinders interacting through the soft spherocylinder model of Lintuvuori and Wilson [23,58]. The
angle between adjacent spherocylinders is set at 150 degrees. A torsional potential about the cen-
tral spherocylinder is used to control the orientation of the two arms to impose an average planar
conformation corresponding to a typical bent-core mesogen. The snapshot shows a NTB phase
that arises spontaneously on cooling from a higher temperature nematic phase. (b): Twisted layer
structures arising from simulations of 10,000 chiral-rigid-asymmetric-bent-core molecules. Green and
white spheres represent a central aromatic core and interact through a Lennard–Jones potential with
σ/σ0 = 0.75, the pink spheres interact through a WCA potential with σ/σ0 = 0.75. The two arms are
positioned at an angle of 120 degrees.
Crystals 2022, 12, 685 8 of 24

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–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, γ1 , which is important (for
example) in determining the on and off times of a twisted nematic display (toff ∝ γ1 ,
ton ∝ γ1 ) and for other nematic devices. γ1 can be determined by a number of techniques
including equilibrium and non-equilibrium molecular dynamics methods [66–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 γ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, γ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 es and eb 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–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 achi-
ral 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 ex-
Crystals 2022, 12, 685 9 of 24

act opposite of what would be expected. Calculations of the chirality order parameter,
χ, show that although these molecules are on average achiral, they can exhibit “chiral
conformations, which have extremely large values of χ in comparison to standard chiral
dopants [83–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π/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, ∆µ, is directly proportional to the helical twisting power, β, and the twist elastic
constant K2 [89].
∆µ = µ− − µ+ = 8πβK2 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 ex-
tended for tens, or (in some cases) hundreds, of nanoseconds. Moreover, coarse-grained
simulation methods have become very useful in studying traditional lyotropic phase dia-
Crystals 2022, 12, 685 10 of 24

grams. 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.

a) b) O
O
S
O
-
R1 +
Na

R2
R1+R2=C11 H24

c)
d)

e)

f)

Figure 6. Structures and simulation models for versions of LAS (linear alkylbenzene sulfonates).
(a) Molecular structure of a fully linear single chained version of the anionic surfactant LAS; (b) chem-
ical structure of typical branched LAS molecules used in industry, (c) single micelle of LAS in water
(sulfonate head groups are shown in yellow and red, and blue sites represent sodium counter ions,
water molecules are shown in a partially transparent representation); (d) three DPD models of LAS
with the orange bead representing the sulfonate head group, the yellow bead representing the phenyl
group, and the purple beads representing parts of the alkyl chain; (e,f) simulations snapshots from
the phase diagram of the linear form of LAS adapted with permission from Ref. [91], 2018, Sarah J.
Gray; showing (e) the hexagonal phase composed of cylindrical micelles and (f) the lamellar phase.
Crystals 2022, 12, 685 11 of 24

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:

FC,ij = aij (1 − rij /rcut ) for|rij | < rcut

(5)
FC,ij = 0 for|rij | ≥ rcut

which leads to a harmonic repulsive potential. The chemical interactions are mainly
controlled by the parameters aij , 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–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 solu-
bilizing 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 ex-
hibit 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 forma-
tion 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.
Crystals 2022, 12, 685 12 of 24

a) b) c)

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) function-
alised 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.

a) b) c)

Figure 8. The nonionic chromonic mesogen TP6EO2M represented by three different levels of models.
(a) All-atom model (b) a MARTINI-style coarse-grained model, and (c) a simple dissipative particle
dynamics model. For coarse-grained models, bonds (not shown) link adjacent sites and angle
interactions help to define molecular shape, blue beads represent different hydrophobic sites, and red
and orange beads represent hydrophilic sites.

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, ∆Gagg◦ , can

be approximated from the depth of the “attractive well” obtained from potential of mean
force (PMF) calculations, UPMF (r ). 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 c , obtained at a
series of separation distances, s:
Z rmax  
2k T
UPMF (r ) = h f c is + B ds (6)
r s
Crystals 2022, 12, 685 13 of 24

where r is the distance for the PMF and 2kB T/s is a kinetic entropy term, which accounts for
the increase in rotational volume at larger separation distances [110–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 distri-
bution of aggregate sizes. For TP6EO2M, the binding energy for a n-mer (∆Gagg ◦ = −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, ∆Gagg ◦ ,

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 ∆Gagg ◦ .

  
◦ /RT
∂ ∆Gagg
∆H ◦ = R . (7)
∂(1/T )

Interestingly, for TP6EO2M, − T∆S◦ and ∆H ◦ both contribute favourably to ∆Gagg ◦

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 ∆Gagg ◦ relative to the DPD energy scale, i.e., e that controls the

magnitude of the reduced temperature kB T/e. 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 simu-
lated for TP6EO2M, producing an approximately exponential distribution of stacks in the
isotropic phase (in-line with quasi-isodesmic association); a nematic N phase at intermedi-
ate 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
Crystals 2022, 12, 685 14 of 24

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–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).

“TOP DOWN” ROUTE

SAFT-γ Mie Martini 3

Hybrid force Iterative

matching Boltzmann
(HFM) Inversion (IBI)

“BOTTOM UP” ROUTE

Figure 9. Routes to an ideal (“Goldilocks”) coarse-grained model for chromonic mesogens using
both bottom-up and top-down methodologies; ideally producing a model that is “just right” for
molecular simulation. A Goldilocks model would capture the correct self-assembly behaviour in
solution in terms of the structure of aggregates and the thermodynamics of self-assembly and would
be sufficiently transferable to be used successfully over a range of concentrations and temperatures
(providing of course these were not too “hot” or too “cold”).

Traditionally, structure-based coarse-graining by methods such as multiscale coarse-

graining (MS-CG) in the form of hybrid force matching (HFM) [120–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 (particu-
larly 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
Crystals 2022, 12, 685 15 of 24

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-γ Mie model [2,125–130], which represents coarse-grained sites by
generalised Mie potentials:
"    λa #
σij λr σij
UMie,ij = Ceij − (8)
r r

where exponents λr and λa govern the ranges of the repulsion and attraction of the potential,
and the following is the case.
q

σii3 σjj3 p
eij = 1 − k ij eii e jj . (9)
σij3

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-γ Mie model for TP6EO2M [114].
The balance of k ij values was found to be crucial. If ethylene oxide units are poorly sol-
vated, 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 chlo-
ride (PIC), pinacyanol chloride (PCYN), 5,50 ,6,60 -tetrachloro-1,10 ,3,30 -tetraethylbenzimidazolyl-
carbocyanine chloride (TTBC), and 1,10 -disulfopropyl-3,30 -diethyl-5,50 ,6,60 -tetrachloro-benz-
imidazolylcarbocyanine 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.
Crystals 2022, 12, 685 16 of 24

Examples of the aggregate structures seen are shown in Figure 11. What is particu-
larly 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–133].

(a) (b)

Figure 10. (a) Snapshot from the N phase composed of 900 TP6EO2M molecules simulated at a
coarse-grained level using the Martini 3 model from reference [113] (54,909 Martini water sites
representing 219,636 water molecules). The snapshot shows short chromonic stacks that are joined
across the periodic boundary conditions. Small blue dots between columns represent Martini water
sites. (b) Snapshot from the M phase of TP6EO2M showing the hexagonal packing of columns in this
phase (9266 Martini water sites).

PIC PCYN TTBC

H J H J
H

Shift

BIC

Figure 11. (Top left): H and J aggregate structures for PIC, plus shift defects and Y-defects seen in
large aggregates. (Top middle): H aggregate structures for PCYN. (Top right): H and J aggregate
structures for TTBC and a brick-like layer structure seen for TTBC at high concentrations, where
J-aggregation dominates over H-aggregation. (Bottom): H and J aggregate structures for BIC and a
layer structure observed for BIC molecules (side and top view). Figure adapted from Ref. [107], 2021
CC-BY licence, with permission from the Royal Society of Chemistry.
Crystals 2022, 12, 685 17 of 24

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, includ-
ing 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-diethyl-
aminoethyl)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.

10
-1
UPMF / kJ mol

-20

-50
0.2 0.8 1.4 2
r / nm

Figure 12. Cont.

Crystals 2022, 12, 685 18 of 24

Figure 12. (Top left): the atomistic structure of PER and a coarse-grained representation suitable for
use with bottom-up coarse-graining methods and the top-down Martini 3 approach. (Top middle):
the potential of mean force for the binding of a PER dimer calculated from atomistic simulations,
together with the lowest energy structure seen at the bottom of the potential well. (Top right): A
snapshot from atomistic simulations showing the typical aggregate structure for PER. (Middle left):
Aggregates of PER from a tuned Martini 3 model. (Middle right): A snapshot from the chromonic
N phase of PER (water molecules not shown) obtained from long simulation runs using the tuned
Martini 3 model developed in reference [8]. (Bottom left): Plan view of the N phase at 30 wt%.
(Bottom right): Plan view of the M phase at 50 wt%. The top row sub-figures are redrawn from
Figures 2 and 3 of reference [8]. Parts of this figure are adapted from reference [8], 2022 CC-BY licence.

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, TNI 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
Crystals 2022, 12, 685 19 of 24

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 Sections 3.3 and 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 consid-
erably. 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.

Author Contributions: Conceptualization, M.R.W.; methodology, M.R.W., G.Y., T.D.P., M.W., S.J.G.,
J.L. and N.J.B.; software, M.R.W., G.Y., T.D.P. and M.W.; validation, M.R.W., G.Y., T.D.P., M.W., S.J.G.,
J.L. and N.J.B.; formal analysis, M.R.W., G.Y., T.D.P., M.W., S.J.G., J.L. and N.J.B.; investigation, M.R.W.,
G.Y., T.D.P., M.W., S.J.G., J.L. and N.J.B.; resources, M.R.W.; data curation, M.R.W., G.Y., T.D.P., M.W.,
S.J.G., J.L. and N.J.B.; writing—original draft preparation, M.R.W.; writing—review and editing,
M.R.W., G.Y., T.D.P., M.W., S.J.G., J.L. and N.J.B.; visualization, M.R.W., G.Y., T.D.P., M.W., S.J.G.,
J.L. and N.J.B.; supervision, M.R.W. and M.W.; project administration, M.R.W.; funding acquisition,
M.R.W. All authors have read and agreed to the published version of the manuscript.
Funding: This research was part-funded by EPSRC in the UK: EPSRC grant EP/R513039/1 funding a
studentship for Gary Yu, grants EP/J004413/1 and EP/P007864/1 providing funding for Dr. Martin
Walker, and a studentship award (grant EP/M507854/1, award reference 1653213) funding Thomas
Potter. Sarah Gray was supported by an EPSRC Case studentship, reference 1482182, with funding
from Procter and Gamble, and Jing Li was partially funded by a research grant from Procter and
Gamble. Prof. Mark Wilson was part funded by EPSRC by grant EP/V056891/1. This work made
use of the facilities of the N8 Centre of Excellence in Computationally Intensive Research (N8 CIR)
provided and funded by the N8 research partnership and EPSRC (Grant No. EP/T022167/1). The
Centre is coordinated by the Universities of Durham, Manchester, and York.
Data Availability Statement: Not applicable.
Acknowledgments: The authors acknowledge support from the Advanced Research Computing
Division at Durham University for the provision of computer time through its Hamilton high
performance computer system. The authors are grateful for computer time on the Bede supercomputer
at the N8 Centre of Excellence in Computationally Intensive Research.
Crystals 2022, 12, 685 20 of 24

Conflicts of Interest: The authors declare no conflicts of interest. The funders had no role in the
design of the study; in the collection, analyses, or interpretation of data; in the writing of the
manuscript; or in the decision to publish the results.

References
1. Güryel, S.; Walker, M.; Geerlings, P.; De Proft, F.; Wilson, M.R. Molecular dynamics simulations of the structure and the
morphology of graphene/polymer nanocomposites. Phys. Chem. Chem. Phys. 2017, 19, 12959–12969. doi:10.1039/C7CP01552F.
2. Tasche, J.; Sabattié, E.F.D.; Thompson, R.L.; Campana, M.; Wilson, M.R. Oligomer/Polymer Blend Phase Diagram and Sur-
face Concentration Profiles for Squalane/Polybutadiene: Experimental Measurements and Predictions from SAFT-γ Mie and
Molecular Dynamics Simulations. Macromolecules 2020, 53, 2299–2309. doi:10.1021/[Link].9b02155.
3. Rodgers, T.L.; Townsend, P.D.; Burnell, D.; Jones, M.L.; Richards, S.A.; McLeish, T.C.B.; Pohl, E.; Wilson, M.R.; Cann, M.J.
Modulation of Global Low-Frequency Motions Underlies Allosteric Regulation: Demonstration in CRP/FNR Family Transcription
Factors. PLoS Biol. 2013, 11, e1001651. doi:10.1371/[Link].1001651.
4. McLeish, T.C.B.; Rodgers, T.L.; Wilson, M.R. Allostery without conformation change: Modelling protein dynamics at multiple
scales. Phys. Biol. 2013, 10, 056004. doi:10.1088/1478-3975/10/5/056004.
5. Marrink, S.J.; de Vries, A.H.; Tieleman, D.P. Lipids on the move: Simulations of membrane pores, domains, stalks and curves.
Biochim. Biophys. Acta (BBA)-Biomembr. 2009, 1788, 149–168. doi:10.1016/[Link].2008.10.006.
6. Catte, A.; White, G.F.; Wilson, M.R.; Oganesyan, V.S. Direct Prediction of EPR Spectra from Lipid Bilayers: Understanding
Structure and Dynamics in Biological Membranes. ChemPhysChem 2018, 19, 2183–2193. doi:10.1002/cphc.201800386.
7. Catte, A.; Wilson, M.R.; Walker, M.; Oganesyan, V.S. Antimicrobial action of the cationic peptide, chrysophsin-3: A coarse-grained
molecular dynamics study. Soft Matter 2018, 14, 2796–2807. doi:10.1039/C7SM02152F.
8. Yu, G.; Wilson, M.R. Molecular simulation studies of self-assembly for a chromonic perylene dye: All-atom studies and new
approaches to coarse-graining. J. Mol. Liq. 2022, 345, 118210. doi:10.1016/[Link].2021.118210.
9. Prasitnok, K.; Wilson, M.R. A coarse-grained model for polyethylene glycol in bulk water and at a water/air interface. Phys.
Chem. Chem. Phys. 2013, 15, 17093–17104. doi:10.1039/C3CP52958D.
10. Anderson, P.M.; Wilson, M.R. Molecular dynamics simulations of amphiphilic graft copolymer molecules at a water/air interface.
J. Chem. Phys. 2004, 121, 8503–8510. doi:10.1063/1.1796251.
11. Allen, M.P. Molecular simulation of liquid crystals. Mol. Phys. 2019, 117, 2391–2417. doi:10.1080/00268976.2019.1612957.
12. Wilson, M.R. Progress in computer simulations of liquid crystals. Int. Rev. Phys. Chem. 2005, 24, 421–455.
13. Wilson, M.R. Molecular simulation of liquid crystals: Progress towards a better understanding of bulk structure and the prediction
of material properties. Chem. Soc. Rev. 2007, 36, 1881–1888. doi:10.1039/b612799c.
14. Care, C.M.; Cleaver, D.J. Computer simulation of liquid crystals. Rep. Prog. Phys. 2005, 68, 2665–2700. doi:10.1088/0034-
4885/68/11/R04.
15. Bremer, M.; Kirsch, P.; Klasen-Memmer, M.; Tarumi, K. The TV in Your Pocket: Development of Liquid-Crystal Materials for the
New Millennium. Angew. Chem. Int. Ed. 2013, 52, 8880–8896. doi:10.1002/anie.201300903.
16. Cheung, D.L.; Clark, S.J.; Wilson, M.R. Calculation of the rotational viscosity of a nematic liquid crystal. Chem. Phys. Lett. 2002,
356, 140–146. doi:10.1016/S0009-2614(02)00380-9.
17. Cheung, D.L.; Clark, S.J.; Wilson, M.R. Calculation of flexoelectric coefficients for a nematic liquid crystal by atomistic simulation.
J. Chem. Phys. 2004, 121, 9131–9139. doi:10.1063/1.1802231.
18. Berardi, R.; Lintuvuori, J.S.; Wilson, M.R.; Zannoni, C. Phase diagram of the uniaxial and biaxial soft–core Gay–Berne model. J.
Chem. Phys. 2011, 135, 134119. doi:10.1063/1.3646310.
19. Allen, M.P.; Warren, M.A.; Wilson, M.R. Molecular-dynamics simulation of the smectic-A* twist grain- boundary phase. Phys.
Rev. E 1998, 57, 5585–5596. doi:10.1103/PhysRevE.57.5585.
20. Lyulin, A.; Al-Barwani, M.; Allen, M.; Wilson, M.; Neelov, I.; Allsopp, N. Molecular dynamics simulation of main chain liquid
crystalline polymers. Macromolecules 1998, 31, 4626–4634. doi:10.1021/ma971105y.
21. Al Sunaidi, A.; Den Otter, W.K.; Clarke, J.H.R. Liquid-crystalline ordering in rod-coil diblock copolymers studied by mesoscale
simulations. Philos. Trans. R. Soc. Lond. A 2004, 362, 1773–1781. doi:10.1098/rsta.2004.1414.
22. Stimson, L.M.; Wilson, M.R. Molecular dynamics simulations of side chain liquid crystal polymer molecules in isotropic and
liquid-crystalline melts. J. Chem. Phys. 2005, 123, 034908. doi:10.1063/1.1948376.
23. Lintuvuori, J.S.; Wilson, M.R. A coarse-grained simulation study of mesophase formation in a series of rod-coil multiblock
copolymers. Phys. Chem. Chem. Phys. 2009, 11, 2116–2125. doi:10.1039/B818616B.
24. Walton, J.; Mottram, N.; McKay, G. Nematic liquid crystal director structures in rectangular regions. Phys. Rev. E 2018, 97, 022702.
doi:10.1103/PhysRevE.97.022702.
25. Al-Barwani, M.S.; Allen, M.P. Isotropic-nematic interface of soft spherocylinders. Phys. Rev. E 2000, 62, 6706–6710.
doi:10.1103/PhysRevE.62.6706.
26. Camp, P.J.; Mason, C.P.; Allen, M.P.; Khare, A.A.; Kofke, D.A. The isotropic-nematic phase transition in uniaxial hard ellipsoid
fluids: Coexistence data and the approach to the Onsager limit. J. Chem. Phys. 1996, 105, 2837–2849. doi:10.1063/1.472146.
27. McGrother, S.C.; Williamson, D.C.; Jackson, G. A re-examination of the phase diagram of hard spherocylinders. J. Chem. Phys.
1996, 104, 6755–6771. doi:10.1063/1.471343.
Crystals 2022, 12, 685 21 of 24

28. Boyd, N.J.; Wilson, M.R. Optimization of the GAFF force field to describe liquid crystal molecules: The path to a dramatic
improvement in transition temperature predictions. Phys. Chem. Chem. Phys. 2015, 17, 24851–24865. doi:10.1039/C5CP03702F.
29. Wang, J.M.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general amber force field. J. Comp.
Chem. 2004, 25, 1157–1174. doi:10.1002/jcc.20035.
30. Onsager, L. The effects of shape on the interaction of colloidal particles. Ann. N. Y. Acad. Sci. 1949, 51, 627–659. doi:10.1111/j.1749-
6632.1949.tb27296.x.
31. Boyd, N.J.; Wilson, M.R. Validating an optimized GAFF force field for liquid crystals: TNI predictions for bent-core meso-
gens and the first atomistic predictions of a dark conglomerate phase. Phys. Chem. Chem. Phys. 2018, 20, 1485–1496.
doi:10.1039/C7CP07496D.
32. Cheung, D.L.; Clark, S.J.; Wilson, M.R. Parametrization and validation of a force field for liquid- crystal forming molecules. Phys.
Rev. E 2002, 65, 051709. doi:10.1103/PhysRevE.65.051709.
33. Wang, J.; Hou, T. Application of Molecular Dynamics Simulations in Molecular Property Prediction. 1. Density and Heat of
Vaporization. J. Chem. Theory Comp. 2011, 7, 2151–2165. doi:10.1021/ct200142z.
34. Caleman, C.; van Maaren, P.J.; Hong, M.; Hub, J.S.; Costa, L.T.; van der Spoel, D. Force Field Benchmark of Organic Liquids:
Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient,
and Dielectric Constant. J. Chem. Theory Comp. 2012, 8, 61–74. doi:10.1021/ct200731v.
35. Qiu, Y.; Smith, D.G.; Boothroyd, S.; Jang, H.; Hahn, D.F.; Wagner, J.; Bannan, C.C.; Gokey, T.; Lim, V.T.; Stern, C.D.; et al.
Development and Benchmarking of Open Force Field v1. 0.0—The Parsley Small-Molecule Force Field. J. Chem. Theory Comput.
2021, 17, 6262–6280. doi:10.1021/[Link].1c00571.
36. Poll, K.; Sims, M.T. An insight into de Vries behaviour of smectic liquid crystals from atomistic molecular dynamics simulations.
J. Mater. Chem. C 2020, 8, 13040–13052. doi:10.1039/D0TC02898C.
37. Poll, K.; Sims, M.T. Sub-layer rationale of anomalous layer-shrinkage from atomistic simulations of a fluorinated mesogen. Mater.
Adv. 2022, 3, 1212–1223.
38. Meyer, R.B. Les Houches Summer School in Theoretical Physics; Balian, R.G., Weil, G., Eds.; Gordon and Breach: New York, NY, USA,
1976; pp. 273–373.
39. Dozov, I. On the spontaneous symmetry breaking in the mesophases of achiral banana-shaped molecules. Europhys. Lett. 2001,
56, 247. doi:10.1209/epl/i2001-00513-x.
40. Chen, D.; Porada, J.H.; Hooper, J.B.; Klittnick, A.; Shen, Y.; Tuchband, M.R.; Korblova, E.; Bedrov, D.; Walba, D.M.; Glaser, M.A.;
Maclennan, J.E.; Clark, N.A. Chiral heliconical ground state of nanoscale pitch in a nematic liquid crystal of achiral molecular
dimers. Proc. Natl. Acad. Sci. USA 2013, 110, 15931–15936. doi:10.1073/pnas.1314654110.
41. Paterson, D.A.; Gao, M.; Kim, Y.K.; Jamali, A.; Finley, K.L.; Robles-Hernández, B.; Diez-Berart, S.; Salud, J.; de la Fuente, M.R.;
Timimi, B.A.; et al. Understanding the twist-bend nematic phase: The characterisation of 1-(4-cyanobiphenyl-40 -yloxy)-6-(4-
cyanobiphenyl-40 -yl)hexane (CB6OCB) and comparison with CB7CB. Soft Matter 2016, 12, 6827–6840. doi:10.1039/C6SM00537C.
42. Cestari, M.; Diez-Berart, S.; Dunmur, D.A.; Ferrarini, A.; de La Fuente, M.R.; Jackson, D.J.B.; Lopez, D.O.; Luckhurst, G.R.; Perez-
Jubindo, M.A.; Richardson, R.M.; et al. Phase behavior and properties of the liquid-crystal dimer 100 ,700 -bis (4-cyanobiphenyl-40 -yl)
heptane: A twist-bend nematic liquid crystal. Phys. Rev. E 2011, 84, 031704. doi:10.1103/PhysRevE.84.031704.
43. Borshch, V.; Kim, Y.K.; Xiang, J.; Gao, M.; Jákli, A.; Panov, V.P.; Vij, J.K.; Imrie, C.T.; Tamba, M.G.; Mehl, G.H.; et al. Nematic
twist-bend phase with nanoscale modulation of molecular orientation. Nat. Commun. 2013, 4, 2635. doi:10.1038/ncomms3635.
44. Yu, G.; Wilson, M.R. All-atom simulations of bent liquid crystal dimers: The twist-bend nematic phase and insights into
conformational chirality. Soft Matter 2022, 18, 3087–3096. doi:10.1039/D2SM00291D.
45. Shadpour, S.; Nemati, A.; Boyd, N.J.; Li, L.; Prévôt, M.E.; Wakerlin, S.L.; Vanegas, J.P.; Salamończyk, M.; Hegmann, E.; Zhu, C.;
et al. Heliconical-layered nanocylinders (HLNCs)–hierarchical self-assembly in a unique B4 phase liquid crystal morphology.
Mater. Horizons 2019, 6, 959–968. doi:10.1039/C9MH00089E.
46. Shadpour, S.; Nemati, A.; Salamończyk, M.; Prévôt, M.E.; Liu, J.; Boyd, N.J.; Wilson, M.R.; Zhu, C.; Hegmann, E.; Jákli, A.I.; et al.
Missing Link between Helical Nano- and Microfilaments in B4 Phase Bent-Core Liquid Crystals, and Deciphering which Chiral
Center Controls the Filament Handedness. Small 2020, 16, 1905591. doi:10.1002/smll.201905591.
47. Berardi, R.; Emerson, A.P.J.; Zannoni, C. Monte Carlo investigations of a Gay—Berne liquid crystal. J. Chem. Soc. Faraday Trans.
1993, 89, 4069–4078. doi:10.1039/FT9938904069.
48. Zannoni, C. Molecular Design and Computer Simulations of Novel Mesophases. J. Mater. Chem. 2001, 11, 2637–2646.
doi:10.1039/B103923G.
49. Berardi, R.; Muccioli, L.; Orlandi, S.; Ricci, M.; Zannoni, C. Computer simulations of biaxial nematics. J. Phys.-Condens. Matter
2008, 20, 463101. doi:10.1088/0953-8984/20/46/463101.
50. Wilson, M.; Allen, M. Computer simulation study of liquid crystal formation in a semi-flexible system of linked hard spheres.
Mol. Phys. 1993, 80, 277–295. doi:10.1080/00268979300102261.
51. Wilson, M.R. Molecular dynamics simulation of semi-flexible mesogens. Mol. Phys. 1994, 81, 675–690. doi:10.1080/00268979400100451.
52. AlSunaidi, A.; den Otter, W.K.; Clarke, J.H.R. Inducement by directional fields of rotational and translational phase ordering in
polymer liquid-crystals. J. Chem. Phys. 2013, 138, 154904. doi:10.1063/1.4798462.
53. Skačej, G.; Zannoni, C. Molecular simulations elucidate electric field actuation in swollen liquid crystal elastomers. Proc. Natl.
Acad. Sci. USA 2012, 109, 10193–10198. doi:10.1073/pnas.1121235109.
Crystals 2022, 12, 685 22 of 24

54. Yan, F.; Hixson, C.A.; Earl, D.J. Self-assembled chiral superstructures composed of rigid achiral molecules and molecular scale
chiral induction by dopants. Phys. Rev. Lett. 2008, 101, 157801. doi:10.1103/PhysRevLett.101.157801.
55. Greco, C.; Ferrarini, A. Entropy-driven chiral order in a system of achiral bent particles. Phys. Rev. Lett. 2015, 115, 147801.
doi:10.1103/PhysRevLett.115.147801.
56. Chiappini, M.; Drwenski, T.; Van Roij, R.; Dijkstra, M. Biaxial, twist-bend, and splay-bend nematic phases of banana-shaped
particles revealed by lifting the “smectic blanket”. Phys. Rev. Lett. 2019, 123, 068001. doi:10.1103/PhysRevLett.123.068001.
57. Fernández-Rico, C.; Chiappini, M.; Yanagishima, T.; de Sousa, H.; Aarts, D.G.; Dijkstra, M.; Dullens, R.P. Shaping colloidal
bananas to reveal biaxial, splay-bend nematic, and smectic phases. Science 2020, 369, 950–955. doi:10.1126/science.abb4536.
58. Lintuvuori, J.S.; Wilson, M.R. A new anisotropic soft-core model for the simulation of liquid crystal mesophases. J. Chem. Phys.
2008, 128, 044906. doi:10.1063/1.2825292.
59. Humpert, A.; Allen, M.P. Elastic constants and dynamics in nematic liquid crystals. Mol. Phys. 2015, 113, 2680–2692.
doi:10.1080/00268976.2015.1067730.
60. Allen, M.P.; Warren, M.A.; Wilson, M.R.; Sauron, A.; Smith, W. Molecular dynamics calculation of elastic constants in Gay-Berne
nematic liquid crystals. J. Chem. Phys. 1996, 105, 2850–2858. doi:10.1063/1.472147.
61. Allen, M.P.; Frenkel, D. Calculation of liquid-crystal Frank constants by computer simulation. Phys. Rev. A 1988, 37, 1813–1816.
doi:10.1103/PhysRevA.37.1813.
62. Fischermeier, E.; Bartuschat, D.; Preclik, T.; Marechal, M.; Mecke, K. Simulation of a hard-spherocylinder liquid crystal with the
pe. Comput. Phys. Commun. 2014, 185, 3156–3161. doi:10.1016/[Link].2014.08.014.
63. Cleaver, D.; Allen, M. Computer Simulations of the elastic properties of liquid crystals. Phys. Rev. A 1991, 43, 1918–1931.
doi:10.1103/PhysRevA.43.1918.
64. Gruhn, T.; Hess, S. Monte Carlo simulation of the director field of a nematic liquid crystal with three elastic coefficients. Z. Naturf.
1996, 51, 1–9. doi:10.1515/zna-1996-1-201.
65. Phuong, N.H.; Germano, G.; Schmid, F. Elastic constants from direct correlation functions in nematic liquid crystals: A computer
simulation study. J. Chem. Phys. 2001, 115, 7227–7234. doi:10.1063/1.1404388.
66. Sarman, S.; Evans, D.J. Statistical mechanics of viscous flow in nematic fluids. J. Chem. Phys. 1993, 99, 9021–9036.
doi:10.1063/1.465570.
67. Kuwajima, S.; Manabe, A. Computing the rotational viscosity of nematic liquid crystals by an atomistic molecular dynamics
simulation. Chem. Phys. Lett. 2000, 332, 105–109. doi:10.1016/S0009-2614(00)01235-5.
68. Sarman, S. Molecular dynamics of liquid crystals. Physica A 1997, 240, 160–172. doi:10.1016/S0378-4371(97)00138-6.
69. Cuetos, A.; Ilnytskyi, J.M.; Wilson, M.R. Rotational viscosities of Gay-Berne mesogens. Mol. Phys. 2002, 100, 3839–3845.
doi:10.1080/0026897021000028410.
70. Osipov, M.; Nemtsov, V. On the statistical theory of the flexoelectric effect in liquid crystals. Sov. Phys. Crystallogr. 1986,
31, 125–130.
71. Ferrarini, A.; Moro, G.; Nordio, P. A shape model for the twisting power of chiral solutes in nematics. Liq. Cryst. 1995, 19, 397–399.
doi:10.1080/02678299508031997.
72. Feltre, L.; Ferrarini, A.; Pacchiele, F.; Nordio, P. Numerical prediction of twisting power for chiral dopants. Mol. Cryst. Liq. Cryst.
Sci. Technol. Sect. A 1996, 290, 109–118. doi:10.1080/10587259608031896.
73. di Matteo, A.; Todd, S.M.; Gottarelli, G.; Solladié, G.; Williams, V.E.; Lemieux, R.P.; Ferrarini, A.; Spada, G.P. Correlation between
molecular structure and helicity of induced chiral nematics in terms of short-range and electrostatic- induction interactions. The
case of chiral biphenyls. J. Am. Chem. Soc. 2001, 123, 7842–7851. doi:10.1021/ja010406r.
74. Ferrarini, A.; Gottarelli, G.; Nordio, P.L.; Spada, G.P. Determination of absolute configuration of helicenes and related biaryls
from calculation of helical twisting powers by the surface chirality model. J. Chem. Soc. Perkin Trans. 2 1999, 411–418.
doi:10.1039/A809593K.
75. Ferrarini, A.; Nordio, P.; Shibaev, P.; Shibaev, V. Twisting power of bridged binaphthol derivatives: Comparison of theory and
experiment. Liq. Cryst. 1998, 24, 219–227. doi:10.1080/026782998207398.
76. Ferrarini, A.; Moro, G.; Nordio, P. Shape model for ordering properties of molecular dopants inducing chiral mesophases. Mol.
Phys. 1996, 87, 485–499. doi:10.1080/00268979600100321.
77. Ferrarini, A.; Moro, G.; Nordio, P. Simple molecular model for induced cholesteric phases. Phys. Rev. E 1996, 53, 681–688.
doi:10.1103/physreve.53.681.
78. Osipov, M.; Pickup, B.; Dunmur, D. A new twist to molecular chirality: Intrinsic chirality indices. Mol. Phys. 1995, 84, 1193–1206.
doi:10.1080/00268979500100831.
79. Solymosi, M.; Low, R.J.; Grayson, M.; Neal, M.P. A generalized scaling of a chiral index for molecules. J. Chem. Phys. 2002,
116, 9875–9881. doi:10.1063/1.1476321.
80. Solymosi, M.; Low, R.J.; Grayson, M.; Neal, M.P.; Wilson, M.R.; Earl, D.J. Scaled chiral indices for ferroelectric liquid crystals.
Ferroelectrics 2002, 277, 483–490. doi:10.1080/00150190214456.
81. Earl, D.J.; Wilson, M.R. Predictions of molecular chirality and helical twisting powers: A theoretical study. J. Chem. Phys. 2003,
119, 10280–10288. doi:10.1063/1.1617980.
82. Neal, M.P.; Solymosi, M.; Wilson, M.R.; Earl, D.J. Helical twisting power and scaled chiral indices. J. Chem. Phys. 2003,
119, 3567–3573. doi:10.1063/1.1590307.
Crystals 2022, 12, 685 23 of 24

83. Earl, D.; Osipov, M.; Takezoe, H.; Takanishi, Y.; Wilson, M. Induced and spontaneous deracemization in bent-core liquid crystal
phases and in other phases doped with bent-core molecules. Phys. Rev. E 2005, 71, 021706. doi:10.1103/PhysRevE.71.021706.
84. Jo, S.Y.; Kim, B.C.; Jeon, S.W.; Bae, J.H.; Walker, M.; Wilson, M.; Choi, S.W.; Takezoe, H. Enhancement of the helical twisting
power with increasing the terminal chain length of nonchiral bent-core molecules doped in a chiral nematic liquid crystal. RSC
Adv. 2017, 7, 1932–1935. doi:10.1039/C6RA27158H.
85. Kim, B.C.; Walker, M.; Jo, S.Y.; Wilson, M.R.; Takezoe, H.; Choi, S.W. Effect of terminal chain length on the helical twisting power
in achiral bent-core molecules doped in a cholesteric liquid crystal. RSC Adv. 2018, 8, 1292–1295. doi:10.1039/C7RA11589J.
86. Lintuvuori, J.S.; Yu, G.; Walker, M.; Wilson, M.R. Emergent chirality in achiral liquid crystals: Insights from molecular simulation
models of the behaviour of bent-core mesogens. Liq. Cryst. 2018, 45, 1996–2009. doi:10.1080/02678292.2018.1492037.
87. Earl, D.J.; Wilson, M.R. Calculations of helical twisting powers from intermolecular torques. J. Chem. Phys. 2004, 120, 9679–9683.
doi:10.1063/1.1718122.
88. Germano, G.; Allen, M.P.; Masters, A.J. Simultaneous calculation of the helical pitch and the twist elastic constant in chiral liquid
crystals from intermolecular torques. J. Chem. Phys. 2002, 116, 9422–9430. doi:10.1063/1.1475747.
89. Allen, M.P. Calculating the helical twisting power of dopants in a liquid crystal by computer simulation. Phys. Rev. E 1993,
47, 4611–4614. doi:10.1103/PhysRevE.47.4611.
90. Wilson, M.R.; Earl, D.J. Calculating the helical twisting power of chiral dopants. J. Mater. Chem. 2001, 11, 2672–2677.
doi:10.1039/B103520G.
91. Gray, S.J. Dissipative Particle Dynamics Simulations of Surfactant Systems: Phase Diagrams, Phases and Self-Assembly.
Ph.D. Thesis, University of Durham, Durham, UK, 2018.
92. Groot, R.D.; Warren, P.B. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem.
Phys. 1997, 107, 4423–4435. doi:10.1063/1.474784.
93. Lavagnini, E.; Cook, J.L.; Warren, P.B.; Hunter, C.A. Translation of Chemical Structure into Dissipative Particle Dynamics
Parameters for Simulation of Surfactant Self-Assembly. J. Phys. Chem. B 2021, 125, 3942–3952. doi:10.1021/[Link].1c00480.
94. Eslami, H.; Khani, M.; Müller-Plathe, F. Gaussian charge distributions for incorporation of electrostatic interactions in
dissipative particle dynamics: Application to self-assembly of surfactants. J. Chem. Theory Comput. 2019, 15, 4197–4207.
doi:10.1021/[Link].9b00174.
95. McDonagh, J.L.; Shkurti, A.; Bray, D.J.; Anderson, R.L.; Pyzer-Knapp, E.O. Utilizing machine learning for efficient parameteriza-
tion of coarse grained molecular force fields. J. Chem. Inf. Model. 2019, 59, 4278–4288. doi:10.1021/[Link].9b00646.
96. Johnston, M.A.; Duff, A.I.; Anderson, R.L.; Swope, W.C. Model for the Simulation of the Cn Em Nonionic Surfactant Family
Derived from Recent Experimental Results. J. Phys. Chem. B 2020, 124, 9701–9721. doi:10.1021/[Link].0c06132.
97. Lydon, J. Chromonic review. J. Mater. Chem. 2010, 20, 10071–10099. doi:10.1039/B926374H.
98. Bosire, R.; Ndaya, D.; Kasi, R.M. Recent progress in functional materials from lyotropic chromonic liquid crystals. Polym. Int.
2021, 70, 938–943. doi:10.1002/pi.6113.
99. Shiyanovskii, S.V.; Lavrentovich, O.D.; Schneider, T.; Ishikawa, T.; Smalyukh, I.I.; Woolverton, C.J.; Niehaus, G.D.; Doane,
K.J. Lyotropic chromonic liquid crystals for biological sensing applications. Mol. Cryst. Liq. Cryst. 2005, 434, 587–598.
doi:10.1080/15421400590957288.
100. Shiyanovskii, S.V.; Schneider, T.; Smalyukh, I.I.; Ishikawa, T.; Niehaus, G.D.; Doane, K.J.; Woolverton, C.J.; Lavrentovich, O.D.
Real-time microbe detection based on director distortions around growing immune complexes in lyotropic chromonic liquid
crystals. Phys. Rev. E 2005, 71, 020702. doi:10.1103/PhysRevE.71.020702.
101. Kaznatcheev, K.V.; Dudin, P.; Lavrentovich, O.D.; Hitchcock, A.P. X-ray microscopy study of chromonic liquid crystal dry film
texture. Phys. Rev. E 2007, 76, 61703. doi:10.1103/PhysRevE.76.061703.
102. Park, H.S.; Agarwal, A.; Kotov, N.A.; Lavrentovich, O.D. Controllable Side-by-Side and End-to-End Assembly of Au Nanorods
by Lyotropic Chromonic Materials. Langmuir 2008, 24, 13833–13837. doi:10.1021/la803363m.
103. Zhou, S.; Sokolov, A.; Lavrentovich, O.D.; Aranson, I.S. Living liquid crystals. Proc. Natl. Acad. Sci. USA 2014, 111, 1265–1270.
doi:10.1073/pnas.1321926111.
104. Walker, M.; Masters, A.J.; Wilson, M.R. Self-assembly and mesophase formation in a non-ionic chromonic liquid crystal system: In-
sights from dissipative particle dynamics simulations. Phys. Chem. Chem. Phys. 2014, 16, 23074–23081. doi:10.1039/C4CP03092C.
105. Walker, M.; Wilson, M.R. Formation of complex self-assembled aggregates in non-ionic chromonics: Dimer and trimer columns,
layer structures and spontaneous chirality. Soft Matter 2016, 12, 8588–8594. doi:10.1039/c6sm01669c.
106. Chami, F.; Wilson, M.R. Molecular Order in a Chromonic Liquid Crystal: A Molecular Simulation Study of the Anionic Azo Dye
Sunset Yellow. J. Am. Chem. Soc. 2010, 132, 7794–7802. doi:10.1021/ja102468g.
107. Yu, G.; Walker, M.; Wilson, M.R. Atomistic simulation studies of ionic cyanine dyes: Self-assembly and aggregate formation in
aqueous solution. Phys. Chem. Chem. Phys. 2021, 23, 6408–6421. doi:10.1039/D0CP06205G.
108. Thind, R.; Walker, M.; Wilson, M.R. Molecular Simulation Studies of Cyanine-Based Chromonic Mesogens: Spontaneous
Symmetry Breaking to Form Chiral Aggregates and the Formation of a Novel Lamellar Structure. Adv. Theory Simul. 2018,
1, 1800088. doi:10.1002/adts.201800088.
109. Carbone, P.; Varzaneh, H.A.K.; Chen, X.; Müller-Plathe, F. Transferability of coarse-grained force fields: The polymer case. J.
Chem. Phys. 2008, 128, 064904. doi:10.1063/1.2829409.
Crystals 2022, 12, 685 24 of 24

110. Villa, A.; Peter, C.; van der Vegt, N.F.A. Self-assembling dipeptides: Conformational sampling in solvent-free coarse-grained
simulation. Phys. Chem. Chem. Phys. 2009, 11, 2077–2086. doi:10.1039/B818144F.
111. Villa, A.; van der Vegt, N.F.A.; Peter, C. Self-assembling dipeptides: Including solvent degrees of freedom in a coarse-grained
model. Phys. Chem. Chem. Phys. 2009, 11, 2068–2076. doi:10.1039/B818146M.
112. Li, C.; Shen, J.; Peter, C.; van der Vegt, N.F.A. A Chemically Accurate Implicit-Solvent Coarse-Grained Model for Polystyrenesul-
fonate Solutions. Macromolecules 2012, 45, 2551–2561. doi:10.1021/ma202490h.
113. Potter, T.D.; Walker, M.; Wilson, M.R. Self-assembly and mesophase formation in a non-ionic chromonic liquid crystal: Insights
from bottom-up and top-down coarse-grained simulation models. Soft Matter 2020, 16, 9488–9498. doi:10.1039/D0SM01157F.
114. Potter, T.D.; Tasche, J.; Barrett, E.L.; Walker, M.; Wilson, M.R. Development of new coarse-grained models for chromonic liquid
crystals: Insights from top-down approaches. Liq. Cryst. 2017, 44, 1979–1989. doi:10.1080/02678292.2017.1342005.
115. Saiz, L.; Klein, M.L. Computer simulation studies of model biological membranes. Acc. Chem. Res. 2002, 35, 482–489.
doi:10.1021/ar010167c.
116. Talandashti, R.; Mehrnejad, F.; Rostamipour, K.; Doustdar, F.; Lavasanifar, A. Molecular Insights into Pore Formation Mechanism,
Membrane Perturbation, and Water Permeation by the Antimicrobial Peptide Pleurocidin: A Combined All-Atom and Coarse-
Grained Molecular Dynamics Simulation Study. J. Phys. Chem. B 2021, 125, 7163–7176. doi:10.1021/[Link].1c01954.
117. Souza, L.M.; Souza, F.R.; Reynaud, F.; Pimentel, A.S. Tuning the hydrophobicity of a coarse grained model of 1, 2-dipalmitoyl-
sn-glycero-3-phosphatidylcholine using the experimental octanol-water partition coefficient. J. Mol. Liq. 2020, 319, 114132.
doi:10.1016/[Link].2020.114132.
118. Bertrand, B.; Garduño-Juárez, R.; Munoz-Garay, C. Estimation of pore dimensions in lipid membranes induced by peptides and
other biomolecules: A review. Biochim. Biophys. Acta Biomembr. 2021, 1863, 183551. doi:10.1016/[Link].2021.183551.
119. Potter, T.D.; Barrett, E.L.; Miller, M.A. Automated Coarse-Grained Mapping Algorithm for the Martini Force Field and Benchmarks
for Membrane–Water Partitioning. J. Chem. Theory Comput. 2021, 17, 5777–5791. doi:10.1021/[Link].1c00322.
120. Potter, T.D.; Tasche, J.; Wilson, M.R. Assessing the transferability of common top-down and bottom-up coarse-grained molecular
models for molecular mixtures. Phys. Chem. Chem. Phys. 2019, 21, 1912–1927. doi:10.1039/c8cp05889j.
121. Izvekov, S.; Voth, G.A. A multiscale coarse-graining method for biomolecular systems. J. Phys. Chem. B 2005, 109, 2469–2473.
doi:10.1021/jp044629q.
122. Noid, W.G.; Chu, J.W.; Ayton, G.S.; Krishna, V.; Izvekov, S.; Voth, G.A.; Das, A.; Andersen, H.C. The multiscale coarse-
graining method. I. A rigorous bridge between atomistic and coarse-grained models. J. Chem. Phys. 2008, 128, 244114.
doi:10.1063/1.2938860.
123. Noid, W.G.; Liu, P.; Wang, Y.; Chu, J.W.; Ayton, G.S.; Izvekov, S.; Andersen, H.C.; Voth, G.A. The multiscale coarse-
graining method. II. Numerical implementation for coarse-grained molecular models. J. Chem. Phys. 2008, 128, 244115.
doi:10.1063/1.2938857.
124. Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. J. Comput. Chem. 2003,
24, 1624–1636. doi:10.1002/jcc.10307.
125. Lafitte, T.; Apostolakou, A.; Avendaño, C.; Galindo, A.; Adjiman, C.S.; Müller, E.A.; Jackson, G. Accurate statistical associating
fluid theory for chain molecules formed from Mie segments. J. Chem. Phys. 2013, 139, 154504. doi:10.1063/1.4819786.
126. Avendaño, C.; Lafitte, T.; Galindo, A.; Adjiman, C.S.; Jackson, G.; Müller, E.A. SAFT-γ Force Field for the Simulation of Molecular
Fluids. 1. A Single-Site Coarse Grained Model of Carbon Dioxide. J. Phys. Chem. B 2011, 115, 11154–11169. doi:10.1021/jp204908d.
127. Avendaño, C.; Lafitte, T.; Adjiman, C.S.; Galindo, A.; Müller, E.A.; Jackson, G. SAFT-γ Force Field for the Simulation of Molecular
Fluids: 2. Coarse-Grained Models of Greenhouse Gases, Refrigerants, and Long Alkanes. J. Phys. Chem. B 2013, 117, 2717–2733.
doi:10.1021/jp306442b.
128. Müller, E.A.; Jackson, G. Force Field Parameters from the SAFT-γ Equation of State for use in Coarse-Grained Molecular
Simulations. Annu. Rev. Chem. Biomol. Eng. 2014, 5, 405–427. doi:10.1146/annurev-chembioeng-061312-103314.
129. Fayaz-Torshizi, M.; Müller, E.A. Coarse-Grained Molecular Simulation of Polymers Supported by the Use of the SAFT-γ Mie
Equation of State. Macromol. Theory Simul. 2021, 31, 2100031. doi:10.1002/mats.202100031.
130. Fayaz-Torshizi, M.; Müller, E.A. Coarse-grained molecular dynamics study of the self-assembly of polyphilic bolaamphiphiles
using the SAFT-γ Mie force field. Mol. Syst. Des. Eng. 2021, 6, 594–608. doi:10.1039/D1ME00021G.
131. von Berlepsch, H.; Böttcher, C.; Dähne, L. Structure of J-Aggregates of Pseudoisocyanine Dye in Aqueous Solution. J. Phys. Chem.
B 2000, 104, 8792–8799. doi:10.1021/jp000085q.
132. Bricker, W.P.; Banal, J.L.; Stone, M.B.; Bathe, M. Molecular model of J-aggregated pseudoisocyanine fibers. J. Chem. Phys. 2018,
149, 024905. doi:10.1063/1.5036656.
133. Kirstein, S.; Daehne, S. J-aggregates of amphiphilic cyanine dyes: Self-organization of artificial light harvesting complexes. Int. J.
Photoenergy 2006, 5, 020363. doi:10.1155/IJP/2006/20363.