Understanding the properties of liquid-crystalline polymers by computational modeling

A topical review of recent theoretical work on the properties of lyotropic solutions and melts containing semiflexible polymers in thermal equilibrium is given, with a focus on the liquid-crystalline and smectic order of these systems in the bulk and under confinement. Starting with a discussion of single chain properties in terms of the Kratky-Porod worm-like chain model and its limitations, extensions along the lines of Onsager’s theory for the isotropic-nematic transition of solutions of hard rods are briefly reviewed. This discussion is followed by a review of recent Molecular Dynamics simulations and classical Density Functional Theory calculations. It is argued that, even in the simplest case of athermal solutions, coarse-grained polymer models must include three lengths: the contour length L, the persistence length ℓp, and the effective chain thickness D. The discussion is then extended to semiflexible polymers in solutions with thermally driven transitions, where the isotropic-isotropic phase separation competes with the isotropic-nematic transition. Basic physical phenomena due to confinement of these systems in thin films with either repulsive or attractive walls are briefly reviewed, and conditions for the formation of strongly ordered surface layers are identified.


Introduction
Semiflexible polymers are of great interest as constituents of materials in a variety of contexts: due to their liquid-crystalline phases they are of basic interest for applications such as displays and other devices [1][2][3][4][5]. While various other materials formed from small anisotropic molecules are also available for these purposes [4,5], polymeric liquid crystals have various practical advantages, such as controlling some of their Frank elastic constants [4][5][6] by suitably chosen molecular weight, and also their production may be cheaper etc. Other applications of liquid-crystalline polymers include the production of fibers with ultrahigh strength [2]. Finally, we state that liquid crystallinity of polymers may also play a role in biological contexts [7,8], but this aspect and other applications are beyond the scope of the present topical review.
We focus here on the understanding of phase diagrams of semiflexible polymers in lyotropic solutions and in corresponding melts, and discuss the properties of such systems both in the bulk and confinement, from the perspective of theory and computer simulation. In principle, the theoretical physics of materials aims to explain the structure-property relationships, starting with the chemical structure [9]. However, for liquid-crystalline polymers, such an approach building a quantitative bridge from the chemistry to macroscopic properties is still in its infancy. Unlike the simulation approaches for liquid-crystalline materials formed from small molecules [10,11], a chemically realistic model development shall not be dealt with here. Disregarding side-chain liquid-crystalline polymers [12] completely, we focus on a coarse-grained description of main-chain liquid-crystalline polymers in terms of mesoscopic parameters, such as their diameter D, persistence length ℓ p , and contour length L. Furthermore, solvent molecules are not included explicitly in the models discussed here, to facilitate the theoretical description and to expedite the simulations. With such an implicit description, the solvent quality can be controlled by the reduced strength  [16].
ε/(k B T) of some effective potential acting between the coarse-grained repeat units, with Boltzmann's constant k B and temperature T. This implicit treatment of solvent effects is possible, since the dynamic response of the considered systems (viscosity, behavior under shear deformation, etc) is completely disregarded here, although we are well aware that these aspects are important for the processing of such materials. Finally, we consider only polymers with linear chemical architecture, disregarding semiflexible ring polymers [13], as well as branched polymers such as bottle-brush polymers, dendronized polymers, etc. To give an overview over the wide variety of naturally occurring and synthetic semiflexible polymers , table 1 provides examples of estimates for ℓ p and D (taken from references [14,15]). While D is typically in the range from 1 to 2 nm (apart from biopolymers such as F-actin), ℓ p varies over a wide range.
In section 2, we shall briefly characterize models and methods, on which the results reviewed in the later sections are based. Section 3 deals with phase diagrams and properties of bulk homogeneous phases of athermal systems, while section 4 discusses a few examples where phase separation into coexisting phases occurs in solvents of varying quality. Section 5 describes inhomogeneities in the structure of liquid-crystalline polymeric systems arising from confinement by planar walls, considering both purely repulsive and attractive surfaces. Finally, section 6 contains a brief summary.

Models and theoretical methods
The properties of semiflexible polymers can be studied using a wide range of numerical techniques, like Monte Carlo (MC) or Molecular Dynamics (MD) simulation methods [17,18] as well as analytic theoretical methods such as Self-Consistent Field Theory [19] and Density Functional Theory (DFT) [20,21]. At the outset we would like to point out that by 'Computational Modeling' in this review we mean numerical approaches for coarse-grained models starting from the molecular scale. Numerical modeling on macroscopic scales, such as the study of topological defects in nematics through numerical methods based on the Landau-de Gennes free energy density [22], lies beyond its scope. Also, simulations addressing phenomena on mesoscopic scales, where molecular information is introduced in a field-theoretic approach [23,24], are out of consideration here. Such hybrid approaches promise great potential for the future, but they have been applied only occasionally in the contexts discussed here.
Any simplified molecular-level model of the considered semiflexible polymers must allow for two length scales that are both large and independent from each other, i.e. the persistence length ℓ p and the contour length L. Figure 1 provides schematic illustrations of some of the models that have been used for the simulation of semiflexible polymers. In Table 2 we have summarized some frequently used symbols and their physical meaning for improving readability.

Lattice models
In early work, for simplicity and computational efficiency, lattice models such as the self-avoiding walk on square and simple cubic lattices (see figure 1(a)) have been very popular (e.g. references ). Chain stiffness was typically introduced through a bond-bending potential that reduces to some extent the number of 90 • kinks on the lattice. The attractive feature of these particularly simple models is that contour lengths up to L = 80000 lattice spacings were accessible for simulations of single chains. However, it must be noted that the lack of small-angle bending deformations implied by such models means that physics, which is important under certain circumstances, is missed. For instance, chains confined by walls [43,44] do not exhibit fluctuations on the scale of the deflection length [46][47][48][49], but can only exhibit propagation along straight lines over distances of the order of ℓ p or larger [43]. Similarly, the lattice model predicts that adsorption of stiff polymers at planar surfaces occurs already for very small adsorption energy of order ε crit ads / (k B T) ∝ ℓ −1 p [26,32,45], while treatments based on the Kratky-Porod (KP) [50] worm-like chain [51][52][53]. While the majority of previous studies focused on the single chain behavior [26,[30][31][32][34][35][36][37][38][39][40][41][42][43][44][45], it is clear that in multi-chain systems even worse lattice artifacts must occur. For example, the director of nematic order is constrained to be oriented along one of the lattice axes. Thus the character of this order is of the Ising type [54] in d = 2 dimensions or of 3-state Potts type [55] in d = 3. Also generalizations of this lattice model where smaller bond angles are admitted, such as the well-known bond fluctuation model [56,57] augmented by a bond-angle potential [37,[58][59][60][61][62], still suffer from the same unphysical description of the nematic phase. Despite these inherent limitations, such lattice models have some attractive features: 'layering' of effective monomers near hard walls is qualitatively described [60][61][62]; large system sizes (such as films bounded by hard walls at distances L z = 200 to 500 lattice spacings apart, with lateral linear dimensions up to 150 lattice spacings) allow a rather accurate estimation of the phase diagram for 'capillary nematization' [63] for such liquid crystal forming polymers [60][61][62].

Tangent hard sphere chains and related models
Alternative approaches that avoid the artifacts of lattice models and are also well-suited for both analytic methods and for MC simulations are based on hard particles in the continuum. Just as a system of hard spheres is a standard model for dense fluids [64,65], small molecule fluids forming nematic and smectic phases can be modeled by systems of hard ellipsoids or hard spherocylinders, etc [66,67]. One can proceed from such models of nematogens to models of liquid-crystalline polymers in several ways, e.g. tangent hard spheres with a bond bending potential [68], or joining hard spherocylinders with an elastic bending potential [69]. In their pioneering work [69] (we shall discuss some of their results in section 3) Dijkstra and Frenkel used spherocylinders of length l sp and diameter D with aspect ratio of l sp /D = 4 and a bending energy where θ sp is the angle between the unit vectors along the axes of subsequent spherocylindrical segments.
Here C is an elastic bending constant which can be chosen very large so that the system approaches the hard rod limit. However, the crossover to the behavior of fully flexible polymers, that could be easily studied with the generalized self-avoiding walk model [38,39], is less conveniently addressed with this model. A two-dimensional variant of this model was studied in the range of persistence lengths ℓ p = C/ (k B T) from 2 to 120 (choosing l sp = 1 as the unit of length) [70], providing evidence for a continuous transition from the isotropic phase to quasi long-range orientational order.

Onsager theory and its extensions
For analytical methods, such as DFT, the tangent hard sphere model is a convenient starting point, and follows closely the classical work of Onsager describing the isotropic-nematic transition in solutions of long Hexagonal order parameter thin hard rods [71] (which is a reasonable model for, e.g. solutions of the tobacco mosaic virus). For homogeneous systems, the molecular density ρ mol (r, ω) is written as a product of the molecular orientational distribution function, f (ω), and the spatial density distribution, ρ mol = N /V. Here, ω is an abbreviation for the two polar angles (θ, ϕ) characterizing the orientation of a molecule, and N is the number of rods (of length L and diameter D) in a system with volume V. Note that this ansatz disregards any non-uniformity of the particle density (unlike DFT of simple liquids [20]), but rather considers only the non-uniformity of molecular orientations. In the isotropic phase, f (ω) is independent of the molecular orientation, and is simply given by f iso = 1/(4π) in d = 3. The Helmholtz free energy per molecule can then be written as the sum of ideal and excess terms: where V exc (ω, ω ′ ) is the excluded volume for two rods with orientations ω and ω ′ , respectively, when only the second virial coefficient is taken into account. Noting that V exc ∝ DL 2 sin γ, where γ is the angle between them, one finds that the excess term, equation (3), scales proportionally to ρ mol L 2 D (or ∝ ρD 2 L, respectively, when we consider a rod as a straight chain of N = L/D tangent spheres of diameter D, at monomer density (2) is of order unity in the isotropic phase, the onset of the nematic order occurs when the excess term becomes of order unity as well, i.e. for ρ i D 3 ∝ D/L. Thus, for long enough rods the transition occurs at a very small density, and hence the approximations made are well justified. Note that the isotropic-nematic transition must be of first order due to general symmetry arguments based on Landau theory [6], which is born out by explicit calculation [71,72]. So in terms of the dimensionless variable c = (π/4)ℓ p D 2 ρ, one finds that the isotropic phase is stable for c < c i = (π/4)ℓ p D 2 ρ i , while the nematic phase is stable for c > c n , and the region c i < c < c n is a two-phase coexistence region. We have replaced here the contour length L by ℓ p (with the understanding that ℓ p = L for strictly rigid rods), because then the variable c defined here can be used for the case of solutions of semiflexible polymer as it stands. Unfortunately, one can not compute either c i , c n or the distribution function in the nematic phase, f n (ω), exactly by analytical methods, but must resort to numerical methods instead [72]. However, a plausible approximation for f n (ω) [73], that is also useful for describing nematic phases of semiflexible polymers, is where θ is the angle between a rod and the 'director' [6] that characterizes the ordering in the nematic phase. The parameter α controls the amount of bond orientational fluctuations in the liquid crystal. Minimizing the total free energy with respect to α yields One can then show that [72] c i ≈ 3.451, c n ≈ 5.122, S n ≈ 0.91, S n being the nematic order parameter at c = c n . Here we recall that orientational order is described by the traceless tensor where the index i labels the rods which are aligned along unit vectors u i . In the nematic phase, the three eigenvalues of Q αβ are Λ 1 = S, Λ 2 = Λ 3 = −S/2, respectively. This approach describing the isotropic-nematic transition for an athermal solution of very thin long hard rods has been generalized by , Odijk [73], and Chen [77] to semiflexible polymers, which were described by the KP worm-like chain model [50]. In the continuum version of this model (see figure 1(c)), a semiflexible polymer is described by a space curve r(s) where s is a coordinate along this curve, from one chain end at s = 0 to the other chain end at s = L. The Hamiltonian of this model is i.e. the potential energy is due to an integral of the local curvature of the polymer along the chain. From equation (8) one can show that tangent unit vectors u(s) along the chain lose their orientation over the persistence length with bond orientational autocorrelation function, C(s). From equation (9) it is straightforward to derive the mean-square distance between monomeric units a distance s apart (note that R(s) =´s The basic idea to combine this description with Onsager's [71] theory is to divide the chain for L ≫ ℓ p into pieces of length ℓ p (in the opposite limit, L ≪ ℓ p , the chain molecule is rod-like anyway). So the analog of the ideal gas term, equation (2), is then an ideal gas of persistent segments, which can reorient independently of each other, F id /(N k B T) ∝ L/ℓ p . The translational part is still ln(ρ mol ), as in equation (2). The excess part then is the product of the density of persistent segments (ρD/ℓ p ) times the number of persistent segments per chain (L/ℓ p ) times the excluded volume between the two persistent segments, which is of order Dℓ 2 p . This rough order of magnitude argument (which is substantiated by explicit calculations [73][74][75][76][77]) yields F exc /(N k B T) ∝ ρLD 2 . The isotropic-nematic transition occurs when F exc is of the same order as F rot id , i.e. for ρ i ∝ (D 2 ℓ p ) −1 . The explicit calculation can then describe the full crossover from the Onsager rod limit to this result in terms of a function f(L/ℓ p ) as [73][74][75][76][77] with f(x ≫ 1) = const and f(x ≪ 1) ∝ 1/x.

Going beyond the 2nd virial coefficient
A key problem of the treatment based on the second virial coefficient is that it is expected to be accurate only for ℓ p /D ≫ 1; otherwise, higher order virial coefficients come into play as well. However, there is no unique way to extend Onsager-style theory [73][74][75][76][77], and one can find many attempts to do this in the literature (e.g. references [14,78,79]). Here we shall focus on the most recent approach [80,81], which finds support from the accompanying MD simulations (we shall describe this methodology briefly below). Egorov et al [80,81] simply rescaled the prefactor of the excluded volume term, so that instead of ρ mol in equation (3) one uses a rescaled density ρ mol a resc . Following Parsons [82] and Lee [83], one bases the choice of a resc on the Carnahan-Starling (CS) equation of state [84] of a hard sphere fluid of N unconnected spherical monomers: where η = ρ(πD 3 /6) is the monomer packing fraction. The DFT results based on equation (12) henceforth will be termed DFT-CS. In the low density limit η → 0, the Onsager-style results [77] are exactly recovered, since a PL resc → 1.
One could argue that the CS equation of state is inappropriate for describing a dense polymeric fluid, since the chain connectivity is completely disregarded in a PL resc . To test for potential issues of this description, Egorov et al [80,81] proposed an alternative rescaling of the prefactor ρ mol in equation (3), based on the Generalized Flory Dimer (GFD) equation of state [85]. The resulting version of the DFT will be termed DFT-GFD. We emphasize that this approximation is expected to be good only for moderately stiff chains, so that we have an isotropic polymeric fluid still as the stable phase at rather high monomer densities. But one finds that this approximation does not reduce to the Onsager form in the low density limit.
A further difficulty concerns the explicit form of V exc (ω, ω ′ ) which is known for hard rods and also for infinitely stiff hard sphere chains (i.e. the 'shish-kebab' model [86]), but not for semiflexible polymers. Accordingly, Egorov et al [80,81] used an empirical expression for V exc (ω, ω ′ ) due to Fynewever and Yethiraj [87] obtained from fitting MC simulation results for two chains interacting with each other, for various choices of chain lengths and stiffnesses. Thus, the semiflexible chains can be regarded as pseudo rods with mutual excluded volume adjusted to account for the bending degrees of freedom of the chains. The order parameter S resulting from these DFT calculations hence characterizes the order of the whole chains, different from the order parameter defined in equation (7) referring to individual bond vectors [88]. In view of these approximations, it is clear that the DFT results of Egorov et al [80,81] for semiflexible polymers are much less rigorous than DFT results available for simple fluids (e.g. references [89,90]). Nevertheless, we include those results here because their usefulness has been established by extensive comparisons with MD simulation results.
A further caveat needs to be made that equations (2)-(12) all refer to d = 3 dimensions only. One can formulate an Onsager-style theory for hard rods in d = 2 dimensions as well, but readily finds that using only the second virial coefficient never becomes self-consistent [91]. In fact, MC studies of corresponding models have revealed Kosterlitz-Thouless [92] type continuous transitions and an algebraic decay of orientational correlation functions rather than true long-range nematic order [93,94]. An interesting aspect needs also to be noted when the KP model, equation (8), is used in d = 2: while in d = 3 there are two transverse directions to the tangent vector u(s), there is only a single one in d = 2. As a consequence, orientational correlations in d = 2 decay slower than in d = 3, and equation (9) needs to be replaced by C(s) = exp[−s/(2ℓ p )] when one keeps the meaning of ℓ p as defined in d = 3. Then in equation (10) ℓ p needs to be replaced by 2ℓ p throughout as well.

Molecular dynamics (MD) simulation methods
One of the most powerful simulation techniques for soft matter systems, and in particular polymers [17,18,65], is based on the numerical integration of the Newtonian equations of motion for the considered many-body systems. While for the tangent hard-sphere model mentioned in section 2.2 above this technique would be difficult to apply, it can be used for the slightly modified 'rattling spheres' model [95]. In that model, the internal molecular motion consists of bonded atoms 'rattling' in their (narrow) confining potentials, with the number and tolerance of bonding constraints controlling the flexibility of the chain. In that first MD simulation of a liquid-crystal forming system of short semiflexible polymers, the three phases -isotropic, nematic, and smectic A-could be identified, at least qualitatively. In more recent work (e.g. references [80,81,[96][97][98][99][100][101][102]) bead-spring models are mostly used. In the case of lyotropic solutions, the effective pair potential between the polymer beads, U MM , is often modeled as a purely repulsive pair potential U R : where r is the distance between two beads, and U LJ (r) is the standard Lennard-Jones (LJ) potential: with interaction strength ε and bead diameter σ, which are often chosen as the MD units of energy and length, respectively. Equation (13) holds both between bonded and non-bonded beads, and the resulting interaction describes polymers dissolved in a good solvent (recall that the solvent molecules are accounted for only implicitly).
In MD simulations, monomers are typically not bonded via rigid constraints, but rather through (harmonic) springs. Many studies [80,81,97,101,102] use the finitely extensible non-linear elastic (FENE) potential [103] with spring constant k = 30 ε/σ 2 and maximum bond length r 0 = 1.5 σ, to prevent unphysical bond crossing [104,105]. This choice for the parameters leads to an effective bond length ℓ b = 0.97 σ ≈ σ between nearest neighbors along the chain. Thus, the model resembles the tangent hard sphere chain model of DFT, but with smooth rather than hard repulsive potentials. When one is interested in smectic phases at large monomer densities [98][99][100], it may be preferable to replace equation (15) by a simple harmonic potential with a large spring constant k ′ ≫ k B T/σ 2 , and an equilibrium bond length of ℓ b < σ. Alternatively, one can also introduce rigid constraints to fix the distance between consecutive beads [98]. In either case, excluded volume interactions between the nearest-and/or next-nearest neighbors along the chain are typically turned off, leading to chains of strongly overlapping spheres with a smoother surface compared to a bead-spring polymer with FENE bonds. For example, in references [99,100] a value of ℓ b = 0.5 σ was employed for U harm , while in reference [98] rigid bonds with ℓ b = 0.6 σ were utilized. When one wishes to consider variable solvent quality, one may use instead of equation (13) a potential of with cutoff radius r c . The parameter T eff controls the effective temperature of the system (or equivalently the quality of the implicit solvent). The advantage of varying T eff over the thermodynamic temperature, T, is that T eff only affects the strength of the attractive monomer-monomer interaction, U A , while leaving the strength of the bond and bending interactions (and thus ℓ p ) unchanged. This model then leads to phase separation into a polymer-rich and a solvent-rich phase, both in the case of fully flexible polymers, and when a bond-bending potential is added to describe stiff chains [106]. A standard choice for this bond-bending potential is [80,81,97,101,102] where ε b controls the chain stiffness, and θ ijk is the angle between two subsequent bond vectors, a i = r j − r i and a j = r k − r j , with r i , r j , and r k indicating the positions of the three subsequent beads of the chain (see figure 1(d)). The bending potential U bend (θ ijk ) is minimized for an angle of θ ijk = 0, which corresponds to the three monomers i, j, and k lying on a straight line. For sufficiently large ε b ≫ k B T, this potential is approximately harmonic in θ ijk , i.e. U bend (θ ijk ) ≈ ε b θ 2 ijk /2, similar to equation (1). However, equation (18) is preferred when one wishes to study the crossover to fully flexible chains where large bond angles (up to almost θ ijk ≈π) can occur. When the potential U MM (r) is only applied to nearest neighbors along the chain, the excluded volume interaction between the beads is essentially shut off, and the resulting model can be viewed as a discretization of the KP model (equation (8)) when ε b is large. Such discrete versions of the KP model have been considered in the context of the adsorption of single semiflexible chains [107], and related somewhat different models have been studied in this context as well (e.g. references [108][109][110]). An estimate of the persistence length, ℓ p , can then be obtained from , the persistence length ℓ p (computed from equation (19)) would be twice as large. As a caveat, we mention that chains adsorbed on a surface still retain some three-dimensional character [107,111]. In that case, ℓ p /ℓ b = κ still holds for nearest neighbor bonds, and only when larger distances s along the chain are considered, a crossover of the effective persistence length, ℓ eff p , from ℓ p to 2ℓ p can be observed with increasing s [107,111], provided that κ ≫ 1 so that excluded volume for the considered distances s does not yet matter.
If not stated otherwise explicitly, the parameters ε and σ of equation (13) have been chosen as the MD units of energy and length, respectively. Further, all simulations have been conducted at T = 1 unless stated otherwise. With monomer (unit) mass, m, we define τ MD = mσ 2 /ε as the intrinsic MD unit of time. As a final note, we emphasize that we consider here cases where the bond length, ℓ b , and the effective thickness of the cylindrical worm-like polymer (D ≈ σ in our case) are of the same order. When one considers worm-like chains with a complicated chemical architecture, such as bottlebrush polymers [112][113][114] or dendronized polymers [115][116][117], also the case D ≈ ℓ p ≫ ℓ b can occur: these polymers then behave like coarse-grained fully flexible polymers on length scales exceeding ℓ p , while equation (9) can apply for s ≪ ℓ p [117]. Thus, a very rich and diverse physical behavior can be captured even in the framework of the simple coarse-grained model description via equations (13)-(18) by appropriate choice of parameters.

Monte Carlo (MC) simulation methods
When one considers the properties of single chains, often specialized Monte Carlo (MC) algorithms are the method of choice, as they allow for highly efficient moves that generate trial configurations within the desired statistical mechanics ensemble [17,18]. MC methods have been routinely employed for reliably estimating the equilibrium properties of very long chains, e.g. the data shown in figure 2(c) was obtained in this way. For systems containing many chains, it depends on the details of the system one wishes to study and on the conditions of interest, whether MC or MD methods are preferable.
MC simulations for semidilute solutions or melts of long semiflexible polymers are possible if the standard MC move, where small random displacements of randomly chosen monomers are attempted, is complemented by the so-called 'slithering snake' move. Here, one randomly chooses an end monomer of a randomly chosen chain, cuts it off, and attempts a transfer to the opposite end of the same chain, where it is added if the Metropolis test [17,18] is successfully passed. To obtain a reasonable acceptance rate, it is advisable to choose the new bond length and bond angle by the force bias algorithm [97]. A similar method has been applied recently to a model with a soft monomeric repulsion [118,119], for studying the elastic constants of a melt of semiflexible polymers.
MC methods have the principle advantage that one can use the grandcanonical ensemble of statistical mechanics, where the chemical potential µ rather than the thermodynamically conjugate variable, i.e. the particle number N , is used as a control variable in the simulation. Alternatively, MC simulations can be carried out in the so-called Gibbs ensemble [120,121], where two boxes, which can exchange both chains and volume, are simulated together to avoid the problem of finite size interfacial effects. For polymeric systems, such MC simulations involve moves where whole chains are inserted into the system or deleted from it. Such techniques are of interest when phase equilibria are studied, such as the phase separation in a dilute or semidilute polymer solution under poor solvent conditions [122][123][124], and also for the isotropic-nematic phase transition of semiflexible polymers under good solvent conditions [68,87]. Unfortunately, the Metropolis acceptance rate for inserting long chains is extremely small [17,18], which limits the practicability of these techniques to only rather short chains, typically in the order of few tens of monomers. Further, no dynamical information can be extracted from traditional MC methods, thus restricting them to the analysis of static equilibrium properties.

Outlook on other models and approaches
At this point we stress that nematic order in solutions and melts of the polymers modeled by equations (13)-(18) arises only from the interplay of polymer stiffness (described by a sufficiently large persistence length ℓ p ) and an isotropic bead-bead interaction, such as equation (13). There are, however, many models where anisotropic pairwise interactions are chosen, which lead to liquid-crystalline order in polymeric systems: e.g. generalizations of the Gay-Berne model [125] for particles with ellipsoidal shape have been used for a model where elongated mesogenic elements are connected to each other by flexible spacers [126]. On a more coarse-grained level (in the spirit of the Maier-Saupe [127] description), one may use a local coupling between neighboring unit vectors u in , u jn ′ of neighboring chains n and n ′ [23,24,128,129]. Many of these studies of semiflexible polymers in thermotropic solution yield valuable insight on the properties of specific polymers, for example poly(3-hexylthiophene) (P3HT) [128]. Such studies are interesting in their own right, but far from the expertise of the present authors. Further, the limited space available for the present topical review has made restrictions on its scope inevitable, and hence all those studies (and related ones) shall not be discussed further. We only note that with such models nematic order may arise already for rather small values of ℓ p /ℓ b ≈ 3 [129], while for models with isotropic repulsive interactions (equation (13)) distinctly stiffer chains are required (e.g. ℓ p /ℓ b ≳ 8).

Single Chains and Dilute Solutions
For stiff polymers under good solvent conditions, the properties of the solution will still reflect the properties of single isolated polymers when the concentration (or the related density ρ of the effective monomeric units, respectively) is small enough so that no coil overlap occurs. To estimate the overlap concentration [130,131], the radii of individual polymer coils must be known. For semiflexible chains with short enough contour lengths L, excluded volume can be neglected and hence (in d = 3 dimensions) equation (10) with s = L then yields the mean-square end-to-end distance describing the crossover from rod-like behavior, R 2 (L) ≈ L 2 , to Gaussian coils, R 2 (L) ≈ 2ℓ p L, where 2 ℓ p can be considered as the 'Kuhn length' of the chain [131]. A similar result holds for the mean-square gyration radius [132]: where n p = L/ℓ p is the number of stiff segments. However, equation (20) is true for an arbitrary large L only for 'phantom chains' without excluded volume. Already, from simple Flory arguments [133,134], one can show that equation (20) breaks down when L exceeds L * = ℓ p n * p = ℓ 3 p /D 2 . Then both R 2 and R 2 g scale as where R * is the radius of a 'blob' containing n * p segments of length ℓ p . Note that equation (21) has been derived using the approximate Flory exponent ν = 3/5 instead of the more accurate value of ν ≈ 0.588 [130,131], but this slight difference will be disregarded here. While the rod to coil crossover simply occurs for n p ≈ 1, the crossover from Gaussian to swollen coils depends on the polymer diameter, This result shows explicitly that for ℓ p ≈ D no regime of Gaussian behavior is expected, as anticipated above. In figure 2, this crossover scaling description is also extended to the single chain structure factor S(q) [45], with wave vector q [130,131]. S(q) is defined as and it is, at least in principle, accessible experimentally via small angle neutron scattering (SANS). Observing the crossover at qℓ p ≈ 1 can be invoked as an experimental method to estimate ℓ p , but the accuracy of this method to estimate ℓ p (as well as other methods based on various dynamic quantities and their theoretical interpretation [135]) is somewhat doubtful, and often the estimates obtained for ℓ p for nominally the same polymer by different methods disagree by up to a factor of 2. An important experimental problem often is polydispersity. On the theoretical side, the computation of S(q) for single semiflexible chains by analytical methods is a formidable problem, even if one restricts the task to the use of the KP model without excluded volume [136]. Thus, figure 2(b) is completely hypothetical, as there are no monodisperse polymers for which S(q) could be accurately measured over 6 decades in q. Even large scale MC simulations for semiflexible self-avoiding walks on a lattice could verify the description contained in figure 2 only by combining data for a broad range of stiffnesses [41]. As reviewed by Hsu et al [41], the analytical description of S(q) for semiflexible polymers is a long-standing and complicated problem. The crossover from Gaussian to swollen chains, showing up in the behavior of the chain radii and in S(q), has also consequences for the decay of the bond orientational autocorrelation function, C(s), equation (9). The exponential decay is expected in the region s < ℓ p , but for s > ℓ p it is replaced by a power law behavior. In the regime s > L * , where excluded volume is fully operative, we have [138] This universal behavior has been verified for the semiflexible self-avoiding walk model [38], with β ≈ 0.824 (β = 4/5 according to Flory theory). In the intermediate regime, ℓ p < s < L * , a gradual crossover from the exponential decay to the power law, equation (23), was observed. The smallness of the 2nd virial coefficient between a pair of segments of length ℓ p and diameter D, which we have already invoked in the context of the Onsager-Khokhlov-Odijk-Chen theory discussed in sections 2.3 and 2.4, shows up also in the osmotic equation of state of the polymer solution [130,131]. This problem was already considered by Khokhlov [139] who predicted that for monomer concentrations ρ above the overlap concentration ρ * and in the good solvent regime the osmotic pressure p osm scales as while for ρ < ρ * one simply has an ideal-gas like behavior, p osm /(k B T) = ρ/N. Scaling considerations [130,131] imply that for ρ = ρ * the density of monomers inside a coil agrees with the overall density, and hence equation (21) yields  [137], as redrawn in reference [41]. In that experimental work, the persistence length was denoted by q and the mean square gyration radius as ⟨S 2 ⟩. Broken straight lines are power law fits with effective exponents, while the full straight line is the Kratky-Porod model. Adapted and reprinted from reference [41], with the permission of AIP Publishing.
For ρ > ρ * it is implied that excluded volume is effective only for length scales less than a screening length (or correlation length, respectively) ξ and the chain can be considered as a random walk of blobs of radius ξ ∝ ℓ 1/5 p D 4/5 g 3/5 , containing g monomers. However, this picture only makes sense if g by far exceeds L * /ℓ p , and hence ξ exceeds R * = ℓ p n * p , cf equation (21). Since the blobs are space filling, we have ρ ∝ g/ξ 3 and thus the condition ξ = R * yields a density ρ * * above which the excluded volume is no longer felt [140], so p osm /(k B T) ∝ ρ 2 and Comparing this result to equation (11), we see that ρ * * is a much smaller monomer density than the density ρ i where nematic order sets in. Finally, we emphasize that the behavior in d = 2 dimensions (applicable for polymers confined in narrow planar slits or strongly adsorbed on planar walls) is completely different from the picture shown in figure 2. There occurs a crossover from the rod regime immediately to the self-avoiding walk regime (with ν = 3/4 in d = 2), so equation (21) is replaced by  [77]. Symbols are the MD results. Reproduced from reference [81] with permission from The Royal Society of Chemistry.
as soon as L exceeds ℓ p [39,134,141,142]. This fact can be understood readily in terms of Flory arguments [134], and has been confirmed both for semiflexible self-avoiding walks and for bead-spring models by extensive simulations [39,141,142].

Nematic and smectic order in concentrated solutions and melts
When the osmotic equation of state is extended to larger densities and pressures, the weakly first order isotropic-nematic transition shows up in the isotherm pressure vs monomer density via a narrow two-phase coexistence region ( figure 3(a)). At still larger densities, the nematic order parameter S can be clearly observed ( figure 3(b)). It was found [80,81] that the MD simulations agree roughly with the DFT predictions, but there was no clear preference for either of the two DFT variants mentioned in section 2.4 (DFT-CS or DFT-GFD, respectively). When the transition occurs at small density ρ, the Onsager-style theory of Chen [77] is in rather good agreement with respect to the prediction of the densities ρ i , ρ n where the ordering occurs; for smaller ε b and/or smaller N, when these transition densities become rather large, this theory [77] becomes quantitatively unreliable, as expected. An important aspect is that the scaled transition densities (ℓ p /D)D 3 ρ i are not simply functions of L/ℓ p only, as proposed in equation (11), but depend on the ratio D/ℓ p as well ( figure 4). It is seen that the dependence is most pronounced when L/ℓ p is large ( figure 4(a)), and for D/ℓ p ≈ 0.010 − 0.025 the predictions of the Onsager-style theory [77] are rather close to the MD results, although for much thicker chains D/ℓ p ≈ 0.10 they fail rather clearly. We also note that comparisons of data and theoretical calculations plotting ρ i ℓ p /D vs ℓ p /L (as done in figure 4(a)) can be found in the literature in various places (e.g. references [14,69]). However, in those earlier works the significance of the parameter D/ℓ p emphasized in references [80,81] had not yet been recognized.
A particular advantage of MD methods is that one can rather straightforwardly carry out simulations in the constant pressure ensemble [65], using rectangular boxes with linear dimensions (L x , L y , L z ) that can fluctuate independently of each other, so that the long range order in the system is not distorted by the use of periodic boundary conditions. In this way, the phase diagram of the model defined in section 2.5 has been studied also for rather large monomer densities, where also smectic and crystalline phases occur (figures 5 and 6). The yellow shaded stripe indicates the two-phase coexistence region predicted by Chen [77], while symbols connected by broken lines are the MD data, and other symbols denote the experimental data for poly(hexylisocyanate) in toluene (PHIC, stars) [143] and poly(yne)platinum in trichloroethylene (PYP, crosses) [144], respectively. Reproduced from reference [81] with permission from The Royal Society of Chemistry.
The smectic phase is characterized by a density modulation in the z-direction (along the nematic director), ρ(z) ∝ cos(2πz/Λ + ϕ). Both the wavelength Λ and the amplitude of this modulation are of interest. For a quantitative characterization, the collective structure factor S coll (q) is considered: In principle, the smectic order parameter, τ , can be defined from the amplitude of the largest peak of S coll (q z ), which is expected to occur at q z = 2π/Λ. However, it was found that global order was often somewhat imperfect due to local defects (e.g. spiral disclinations). Thus, one could obtain better results by computing τ in rectangular slabs with volume l × l × L z , with l = 5 or l = 10, instead of the entire system. The order parameter τ obtained in this way is shown in figure 6(c), as well as an order parameter Ψ 6 defined from which characterizes the transverse order in the (xy)-planes perpendicular to the director. Labeling the chains with index α in an individual smectic layer (M being the number of chains in the layer), one asks whether the center of mass coordinates r CM of the chains projected onto the xy-plane form a triangular lattice. In equation (29), ϕ αβ is the angle between the vector r CM ⊥,β − r CM ⊥,α and the x-axis r CM α ≡ (r CM ⊥,α , z CM α ). A nonzero order parameter ⟨Ψ 6 ⟩ hence indicates that the stretched out rod-like chains form a triangular lattice in their perpendicular plane. Indeed, figure 6 provides clear evidence for a continuous nematic-smectic transition of the considered model, while the transition from the smectic phase to the crystal is clearly of first order. As a caveat, we mention that in real polymeric systems the smectic phase is easily suppressed by the polydispersity of the chains. For a description of possible crystal phases of polymeric systems, a chemically detailed description of the dense packing of the atoms will be required; while we maintain that our simplified coarse-grained model should provide a qualitatively reasonable description of nematic order, no such claim is made concerning the smectic and crystal phases at high density.
A particular advantage of MD simulations of nematic phases of semiflexible macromolecules is that one can elucidate their structure in great detail. For instance, one finds that the persistence length ℓ p , extracted from the angle θ ijk between subsequent bond vectors via equation (19), does not depend on the monomer density of the solution (and hence the local environment encountered by the bond vector), as long as one stays in the isotropic phase. However, the transition to the nematic phase leads to a (discontinuous) increase of ℓ p , and the deeper one moves into the nematic phase the more the progressive alignment of bonds along the director also has an effect of decreasing θ 2 ijk , figure 7(a). One can also investigate the extent to which hairpin conformations are present ( figure 7(b)) and test the corresponding theoretical predictions [145], which were found to be qualitatively reasonable. Another quantity of interest is the so-called 'deflection length' λ defined as λ = ℓ p /α [145]. For polymers with L ≫ ℓ p this length describes the undulations of the local chain orientation around the director (see schematic in figure 8(a)). Figure 8(b) shows a comparison between MD and DFT results for this length, demonstrating excellent agreement between the two methods.  Also the dependence of the parameter α (equation (4)) which dictates the amount of bond orientational fluctuations in the nematic phase is of interest (figure 9(a)). It is seen that the corresponding theoretical prediction has the right order of magnitude, but increases with increasing density much less steeply than the simulation results. Since equation (19) implies θ 2 = 2/α and for large α hence θ 2 is very small, S = (3 cos 2 θ − 1)/2 ≈ 1 − 3 2 θ 2 . When hairpins are neglected, the end-to-end distance is simply given by ⟨R 2 ⟩ ≈ L(1 − θ 2 /2), and one can conclude that 1 − S = 3(1 − ⟨R 2 ⟩/L) in this limit. Figure 9(b) tests for deviations from this simple limiting case, and one can see that the data for ℓ p /L ≈ 1/2 are already close to this limit. On the other hand, for ℓ p /L > 1 the data tend over to a straight line parallel to the ordinate (the latter would be reached in the hard rod limit ℓ p /L → ∞).
We can estimate the extent to which chains are spatially correlated, by defining the radius r ρ from the condition that a (fictitious) cylindrical tube with a chain in it does not contain monomers of any other chains (see figure 8(a)). Since then the monomer density in the tube must coincide with the overall density ρ, ρ = N/(Lπr 2 ρ ) (where Lπr 2 ρ is the tube volume), one concludes that r ρ = (N/(πρL)) 1/2 . For ℓ p /L > 1, the typical angles ⟨θ 2 ⟩ indicate collective fluctuations of neighboring chains that are much larger than expected from r ρ . The collective orientational fluctuations are assumed to be responsible for the fact that the order parameter S saturates with increasing density, according to the simulations, much slower than predicted by the DFT calculations ( figure 3(b)).
In this section, we have only addressed 'monodomain' samples, where the nematic or smectic order is uniform, i.e. free of defects. However, defects are extremely widespread and common in experimental realizations of liquid crystalline systems (see e.g. references [146][147][148][149][150]). Such defects can form when the  system did not have enough time during its preparation to develop perfect order, which can also occur in simulations as a result of improper equilibration or inappropriate starting configurations [102]. Alternatively, defects can emerge due to the influence of external fields, e.g. electromagnetic fields, mechanical stress, and/or external walls. Particularly interesting from an application point of view are walls, which can be exploited to control the local orientation of bonds near the walls through, e.g. anchoring of chains at functional surfaces. Simulations have the advantage that such external factors can be precisely controlled. For example, external walls can be easily incorporated into MD simulations or DFT calculations (see section 5), or they can be completely eliminated through the use of periodic boundary conditions. Indeed, defects and defect-related phenomena play an important role for liquid crystalline materials, but they are out of scope of the present treatment, which focuses on the intrinsic bulk properties of (perfectly) ordered systems.

The Frank elastic constants
In the nematically ordered phase, the resistance of the system to deformations of this nematic order on large length scales is described by the Frank elastic constants K 1 , K 2 , and K 3 [4,6]. Each of these three constants represents a type of distortion of a nematic. The first represents pure splay, the second pure twist, and the third pure bend (see figure 10), and a combination of these constants can be used to represent an arbitrary deformation in a liquid crystal.
Understanding the behavior of these constants is crucial for many applications of liquid-crystalline materials, but it is rather difficult to make reliable predictions of them from either theory or simulations. The free energy density of the nematic, f(r), can be described exclusively in terms of the director field n(r), Equation (30) assumes that, the nematic order parameter S, defined as the largest eigenvalue of the tensor Q αβ (see equation (7)), is close to unity. This is, however, often not the case for liquid crystalline polymers, and also not in our simulations (cf figure 3(b)); for such cases where the nematic order is not perfect, it is better to base the theory on the full tensorial local order parameter [146,147], as considered below in equation (31). Such a full description is indispensable when defects in liquid crystal structures need consideration, which is frequently the case in experiments containing polydomain structures [146]. For the sake of simplicity, we will not give full details on the resulting theoretical description of f(r), but focus on defect-free systems with strong nematic order. One proceeds by assuming that the local orientation n(r) of the director deviates from the mean orientation n 0 only by a small amount, n(r) = n 0 + δn(r), where the orientation of n 0 defines the z-axis. One then defines the Fourier transform of the tensorial order parameterQ αβ (r) (cf equation (7)) as Here, V is the volume of the considered system, and the sums extend over all unit vectors from monomer at site r n,i to the next monomer in the same chain (labeled by index n) at site r n,i+1 . Equation (30) can then be used to compute Q αβ (q) 2 which yields with K R 1 being a wavevector-dependent generalization of K 1 which will not be discussed further here. Considering special orientations of q, e.g. q = (0, q y , q z ), one can separate twist-bend fluctuations from splay-bend fluctuations [101,129]. While equation (32) seems fairly complicated, it must be realized that it is strongly simplified since only orientational fluctuations are considered while packing density fluctuations are completely disregarded. However, the latter must be included when one wishes to address the collective structure factor S(q) as it is accessible in the standard small-angle scattering of neutrons and X-rays. Kamien et al have shown that the structure factor has the following form in the small wave vector limit [151] with mean areal chain packing density ρ 0 , two-dimensional bulk modulus B, and elastic modulus G = k B TL/(2ρ 0 ). This result has in fact been useful for both experimental X-ray scattering studies from nematic polymers [152] and for the interpretation of simulations [97,129].
A key question is to understand how the Frank constants K 1 , K 2 , and K 3 depend on the chain parameters L and ℓ p and the monomer density ρ. This problem is still rather controversially discussed in the literature (see reference [153] and the discussion in reference [101]). In particular, many theories (e.g. reference [154]) imply K 1 = 3K 2 , but this result is derived when one takes orientational fluctuations of the bonds as the only relevant degree of freedom, disregarding the coupling between local density fluctuations and orientational fluctuations completely. Considering the limit L ≫ ℓ p , Vroege and Odijk presented various arguments to suggest that [145] where α controls the amount of bond orientational fluctuations in the nematic phase (see equation (4)). The DFT approach presented by Milchev et al [101] uses again the excluded volume interactions computed by Fynewever and Yethiraj [87] and combines them with an approach due to Lee [155], based on the rescaling of the Onsager-like direct correlation function, similar in spirit to the rescaling of the free energy [83]. This approach had been originally developed for a system of hard spherocylinders, and also  suffers from the problem that it predicts K 1 = 3K 2 . The numerical results thus obtained exhibit a rapid increase of both K 2 and K 3 with density ρ and chain length N ( figure 11(a) and (b)).
The MD results, however, suggest a less dramatic variation with density and chain length (see figure 12). Those results for the elastic constants are in qualitative accord with the results proposed by Vroege and Odijk [145], in particular when one takes into account that α has not reached its limiting value for large N yet for the chosen parameters. Thus, the theoretical situation concerning the behavior of the elastic constants is still somewhat unsatisfactory. Figure 13 shows experimental results for poly-γ-benzylglutamate (PBG) [156,157]. The experiments imply that the twist elastic constant K 2 depends on neither the chain length nor the density (or, equivalently, the volume fraction of the polymers) [156,157]. Indeed, the MD data ( figure 12(b)) show that the chain length dependence of K 2 is much weaker than that of either K 1 or K 3 (note that the ordinate of figure 12(b) extends over a single decade, whereas K 1 and K 3 vary over two and three decades, respectively). In contrast,  [157]. Data extracted by quasielastic Rayleigh scattering. PBG chains dissolved in a mixture of dioxane, dichloromethane, and dimethylformamide to suppress cholesteric helicity and aggregation of PBG chains. Both K1 and K3 were estimated assuming K2 = 0.6 × 10 −7 dyn, independent of ϕ and L/D. Volume fractions in (a)-(c) were deduced from measurements of optical birefringence [156], assuming a linear dependence of the specific birefringence on the volume fraction. Panel (a) reprinted with permission from reference [156]. Panels (b, c) reprinted with permission from reference [157]. Copyright (1988) by the American Physical Society. the DFT results imply a strong chain length dependence of K 2 (cf figure 11(a)), and also the density dependence is much more pronounced. The theory by Vroege and Odijk also implies a non-negligible dependence of K 2 on ρ [145], due to √ α ∝ ρ 1/3 . Unfortunately, the MD data [101] are only available over such a restricted density range that neither the linear density dependence of K 1 ( figure 13(a)), which is supported by equation (34), nor the claim that K 2 is independent of density could be tested so far.
Experiments on the Frank elastic constants of liquid-crystalline polymers are extremely scarce (the methods applied in references [156,157] require the preparation of 'monodomain' samples, i.e. uniform nematic order throughout the sample). Unfortunately, liquid crystalline polymers typically develop non-uniform orientational order, due to disclinations and defects, leading to polydomain structures, rendering the application of formulas such as equations (32,33) impossible. This difficulty explains why experimental results on elastic constants of liquid crystalline polymers are so rare. This lack of experimental data is rather regrettable, since the elastic anisotropy of liquid-crystalline polymers can be much larger than for their small molecule counterparts [156], which could be an advantage with respect to applications. Very large elastic anisotropies have been also found for thermotropic liquid-crystalline polymers by a step-rotation rheo-nuclear magnetic resonance experiment [158]. These anisotropies were observed for a single molecular weight and using melt densities, finding that K 1 exceeds K 2 ≈ K 3 by two orders of magnitude (while the latter are already two orders of magnitude larger than those of low molecular weight liquid crystals).
An interesting lyotropic material is filamentous actin (F-actin), with diameter D ≈ 8 nm and persistence length ℓ p in the range from 10 µm to 17 µm [160,161]. The contour length L can be tuned (in a range comparable to ℓ p ) by suitable capping proteins [162]. This material is of interest for biological or medical applications since it is non-toxic [5]. Zhang et al [163] developed a method to extract ratios of the elastic constants by optical observations of topological defects of thin films of suspensions of F-actin. They analyzed the resulting structures within Landau-de Gennes theory [6], and showed that the ratio K 3 /K 1 increases from about 0.51 to about 2.13, as the contour length increases from L = 1 µm to 2 µm. However, this approach is beyond the scope of the present review, and thus we refer the reader to reference [163] for further details.
Insight on the elastic anisotropy can also be gained from the analysis of the small angle scattering intensity based on equation (33). Note that for isotropic systems S(q)/(k B T) ≈ ρχ T /(1 + q 2 ξ 2 ), where χ T is the isothermal compressibility and ξ plays the role of a correlation length of local density fluctuations. So contours of constant S(q) in the (q z , q ⊥ ) plane are just circles around the origin. Equation (33), however, predicts rather anisotropic contours of constant S(q z , q ⊥ ), see figure 14 [101,152]. Note that the peaks of scattering intensity are offset from the origin. However, neither the experiment [152] nor the simulation [101] did attempt to obtain elastic constants from fitting equation (33) to these noisy data. Clearly, much more effort would be required for making progress along such lines.

Phase separation
Already for flexible chains, one can observe phase separation between a dilute solution and a concentrated solution for temperatures below the critical temperature [106,[164][165][166]. This type of phase separation for isotropic systems of flexible polymers has been extensively studied before, and is not of interest here. Phase separation in solutions containing semiflexible polymers may be driven by entropy or by attractive effective interactions between the monomers (i.e. rather than good solvent conditions one deals with solutions at temperatures below the theta temperature [130,131]). Entropically driven phase separation in systems with purely repulsive effective interactions among the monomers can be expected when we consider, e.g. binary mixtures of semiflexible polymers that differ strongly either in their chain lengths (N 1 ≫ N 2 ) or in their stiffness (ℓ p,1 ≫ ℓ p,2 ) or both [122]. In principle, both isotropic-isotropic and nematic-nematic phase separations are possible, competing with the isotropic-nematic phase separation. While interesting phase diagrams have been predicted by theory [122], this problem has not yet been extensively studied by molecular simulations and hence will not be considered here further. We shall also not discuss (thermotropic) solutions, where the interactions between monomeric units depend on the orientations of the corresponding bond vectors.
In this review, we focus on the case of monodisperse solutions of semiflexible polymers in solvents of varying quality [106,123,124,167], where the monomer-monomer interaction has an additional attractive contribution. Phase separation in these systems has been studied both for lattice models [124] and for bead-spring type models in the continuum [106,123,167]. Padmanabhan et al [167] used MC methods for a model where attractive interactions were of square-well type. Isotropic-isotropic (I-I) phase separation was observed for polymers with chain lengths N = 5 to 30, while for 10 ≤ N ≤ 35 also isotropic-nematic (I-N) transitions could be located. The phase diagrams then have a 'chimney'-type topology. While the isotropic-isotropic critical phenomena are believed to remain in the Ising model universality class (like liquid-vapor phase separation for systems of small molecules [54]), a nontrivial question concerns the stiffness dependence of the (effective) critical temperature T cr eff : while Sheng et al found that stiffness leads to a decrease of T cr eff in comparison with the corresponding flexible chain model [123], other studies reached the opposite conclusion [106,124]. It must be noted, however, that unlike all other studies [106,124], Sheng et al [123] used also a torsional potential in their microscopic model of the semiflexible chains. Thus details of the chosen models are important, and hence one must be careful with conclusions about real systems, where chemical details may matter.
The most recent simulation study by Midya et al [106] presented a detailed comparison between MD and DFT results. In that work, the attraction strength between monomeric units was controlled by the (effective) inverse temperature T −1 eff (see section 2.5). Note that k B T = ε was kept in the MD simulations, so that the spring constant between neighboring beads along the chain described by equations (15) and (18) remained precisely the same as in the purely repulsive case. Another advantage of this model is that it allows for a gradual transition to good solvent conditions by choosing T eff → ∞. In the DFT calculations, the temperature T was varied while simultaneously adjusting the stiffness parameter ε b to achieve a constant reduced stiffness parameter κ ≡ ε b /(k B T).
In reference [106], both the variation of chain length (8 ≤ N ≤ 32) and stiffness (from flexible chains, κ = 0, to semiflexible ones κ = 32) was studied. Figure 15 gives an example of results. One can see that the critical point of I-I phase separation occurs in this model at much lower monomer densities ρ compared to the I-N coexistence, for the chosen parameters. At the temperature of the triple point, the width of the I-N two-phase coexistence region is distinctly wider than at high temperatures. When N and/or κ increases, this width decreases, however. Note that when N → ∞, we expect that the critical density ρ cr → 0 and that the critical temperature T cr eff saturates at the theta temperature [130,131]. On the other hand, when κ increases at fixed N, we expect that ρ cr (κ → ∞) = ρ cr ∞ > 0. From the DFT results ( figure 15(d)) we can expect that the data for ρ cr (κ = 32) and T cr (κ = 32) are already rather close to their limiting values for hard rods (κ = ∞).
The situation is rather intricate, since for increasing chain stiffness the topology of the phase diagram changes: triple point and critical temperature become identical at some specific value κ * (which for N = 32 is slightly larger than κ = 64, according to DFT). For κ > κ * there is a single first order transition from isotropic to nematic solution, and the phase diagram then has a 'swan-neck' topology ( figure 16(b)). The phase separation between the vapor-like (low density) phase and the high density isotropic liquid then occurs as a metastable equilibrium within the vapor-nematic (V-N) two-phase coexistence region only ( figure 16(b)). Due to its mean-field character, DFT yields well-defined metastable states, which in MD simulations (as well as in real systems) have at best a limited 'lifetime' and their observability is doubtful. Another consequence of this mean-field character is the prediction that the shape of the vapor-isotropic (V-I) coexistence curve is parabolic, i.e. ρ I − ρ V ∝ (1 − T/T cr ) β with β = 1/2, while the MD results are compatible with the expected Ising-like behavior [106], β ≈ 0.325. However, the advantage of DFT is that the statistical mechanics is also accessible in ensembles that use intensive variables only, see e.g. figure 16(d), that give a particularly clear description of the possible phase transitions.

Semiflexible chains in confinement
While confinement effects for fully flexible polymers are predominantly entropic [130], confining semiflexible chains incurs both energetic and entropic penalties, where the balance between the two contributions is primarily controlled by the stiffness of the chain and the length scale of the confinement [168,169]. In a nematic system under confinement, the natural tendency of forming an ordered bulk configuration is impeded by the boundary conditions or curvature of the enclosing walls. There are numerous parameters that control the packaging of semiflexible polymers in confinement, such as the size, shape and topology of the confining geometry, the architecture and size of the enclosed polymers, etc. For example, the length scale regime where the contour length of the confined particles is much smaller than the (spherical) container is typical for self-assembled membranes [170] and nematic droplets [169][170][171][172][173][174]. In contrast, the opposite extreme is of particular interest for a range of biological systems, such as DNA-histone complexation in the eukaryotic nucleus [175][176][177] or the packaging of DNA in bacteriophage [178][179][180][181][182][183][184]. Recently, cases where the contour and persistence lengths of the polymers are comparable to the size of the container have been explored as well [185][186][187][188][189][190], revealing interesting phenomena like confinement-induced smectic ordering and topological defects on the container surface. From this short list of examples it is clear that a full review covering previous work on semiflexible polymers in confinement is beyond the scope of this article, and we will focus here instead on the most elementary case of semiflexible chains in planar confinement.
When semiflexible polymers are confined in a slit pore with impenetrable repulsive walls placed with their normals along the x axis, then the polymers experience an 'ordering field' near the confining surfaces that enhances orientations of the bonds parallel to the surfaces. As a consequence, both the nematic order parameter, S(x), and the monomer density, ρ(x), depend in a nontrivial way on the coordinate x perpendicular to the surfaces. Systems in planar confinement were first studied by Escobedo and de Pablo using MC simulations of athermal hard rods [68], where they observed capillary nematization of the confined particles. Large scale simulations using the bond fluctuation model were presented by Ivanov et al [60][61][62], attempting to clarify the confinement-induced order for semiflexible polymers. However, since that model allows orientations of the director only along the lattice axes, one has to be cautious in extending conclusions from those works to real systems.
More recently, off-lattice DFT calculations and MD simulations were performed by some of us [191][192][193][194][195], using the models described in section 2.4 and section 2.5, respectively. Those works considered perfectly flat walls without any structure in the tangential directions, unlike atomistic models where, e.g. the corrugation due to the periodic arrangement of carbon atoms in confining graphite surfaces is taken into account [196].
In the DFT treatment, repulsive walls are simply represented by hard walls located at x = 0 and x = L x , while in the MD framework an additional external potential is used where x ′ denotes the normal distance from the closest confining wall. For simplicity, ε w = ε and σ w = σ was used throughout. Attractive walls were modeled using the Mie 10-4 potential The attraction strength of U Mie was typically set to ε Mie = 3 in references [194,195], and the potential was cut at x ′ = 5. A potential of the type U Mie is appropriate as a coarse-grained description of a surface where the attractive interaction results from a two-dimensional surface layer, e.g. suitable surfactants at a silicon surface can control the wettability of a surface over the full range from complete wetting to drying [197]. . Reproduced from reference [194] with permission from The Royal Society of Chemistry.
In the MD simulations, typically the case of one wall being attractive and the opposite wall being repulsive was simulated. If both walls were attractive, it may happen that the number of chains that adsorb from the solution on the wall at x = 0 may differ from the number of chains adsorbed at the wall at x = L x , although on average they should be equal. Note that one can not simply symmetrize a posteriori the simulated density profiles by mirroring the system along the origin, as some properties such as the nematic order parameter, depend on the local density in a non-linear fashion.
DFT calculations can be carried out rather easily in the grand-canonical ensemble, so that for large enough slit widths L x , the local density in the center of the slit is identical to the density of the corresponding bulk system, i.e. ρ(L x /2) = ρ b . This is, however, not necessarily true for MD simulations in the canonical ensemble, where the total density is strictly conserved, and hence where ρ att and ρ rep describe the excess densities at both walls: This problem of finite size effects of order 1/L x must be carefully considered when one compares results from MD simulations with the corresponding results from DFT calculations ( figure 17). An important quantity for confined systems is the excess free energy due to the walls, i.e. the surface tension. In the MD framework, it can be extracted from the anisotropy of the pressure tensor. In DFT, it follows as an 1/L x correction to the grand potential Ω that is computed [198]. (In the DFT framework, instead of working with a slit of width L x , one can also compute surface tension by considering a semi-infinite system with a single wall and by imposing bulk boundary condition sufficiently far away from the wall.) Figure 18 shows typical results: in the purely repulsive case, the surface tension originates from a subtle balance of various entropic contributions and it is numerically small; in the attractive case, it is orders of magnitude larger and strongly negative, because the chains aligning at the wall lower the potential energy of the system.
When two repulsive walls are chosen, the difference between ρ and ρ b is rather small, and hence can be ignored. Figure 19(a) shows the nematic order parameter S for various choices of L x , and it is evident that the walls stabilize some order already in the isotropic phase. This picture is corroborated by DFT calculations ( figure 19(b)), and hence one can conclude that there is capillary nematization, i.e. the transition to uniform nematic order in the center of the slit occurs at a smaller density ρ b than in the bulk. When the distance between the walls then gets small enough, the discontinuous isotropic-nematic transition disappears in favor of a gradual, non-singular, onset of nematic order with increasing density. This problem has already been studied with a lattice model [60][61][62].
Finally we note that for the strongly attractive wall potential, equation (38), the enhancement of the local density near the wall (figure 17) induces some quasi-two-dimensional order in the adsorbed surface layer, as shown in figure 20 for semiflexible chains with N = 32 [194,195]. While in the case of κ = 3 one finds a standard isotropic fluid ( figure 20(a)), a lot of small nematic clusters form at κ b = 5, including a significant Figure 18. (a) Surface tension as a function of chain length N at bulk density ρ = 0.062 5 and bending stiffness κ = 16 (upper part) and κ = 100 (lower part). Dots are MD results, full curves are from DFT. In the lower part, the DFT result is decomposed into the isotropic part (dashed curve) and the orientational contribution (dash-dotted curve). Adapted and reprinted from reference [191], with the permission of AIP Publishing. (b) Surface tension γ iw of semiflexible polymers in the isotropic phase with N = 32, ρ = 0.1, and the wall potential equation (38) with ε Mie = 3 plotted vs κ. Reproduced from reference [194] with permission from The Royal Society of Chemistry. Figure 19. (a) Variation of the average nematic order parameter (from MD) in the slit with monomer density ρ for N = 32 and κ = 32, for several choices of Lx as indicated. (b) Same as panel (a), but from DFT, and choosing the chemical potential rather than the monomer density as control variable. Reproduced from reference [193] with permission from John Wiley and Sons. fraction of hairpin-like conformations ( figure 20(b)). As the stiffness is further increased to κ = 15, the (defective) nematic order extends throughout the system, as shown in figure 20(c): note, however, that nematic long-range order is ruled out in strictly two-dimensional systems. For the highest stiffness κ = 128, some local order of smectic C type clearly seems to be present ( figure 20(d)). Also strictly long-range smectic order in two dimensions is not expected either, nevertheless systems such as shown in figure 20(d) could find potential applications to manufacture patterned substrates.

Summary
In this topical review, computer simulations and Density Functional Theory (DFT) treatments of semiflexible polymers that can form liquid crystalline structures were discussed. The considered coarse-grained models did not address the specific aspects related to the chemical composition of particular materials, but rather aimed at providing qualitative insights into the generic features of a broad class of materials. Thus, the persistence lengths of the examples mentioned in this review range from a few nm up to 17 µm (in the case of F-actin), depending on the applied mapping.  [194] with permission from The Royal Society of Chemistry.
We have emphasized that a minimal model of semiflexible polymers in this context must be described in terms of three lengths, i.e. the persistence length ℓ p that describes the local bending stiffness, the effective diameter D of the 'worm-like' chain, and the contour length L. For bead-spring models, a fourth length is given by the distance ℓ b between effective monomers along the chain. These subunits are typically much larger than actual 'chemical monomers' , since often it is convenient to choose ℓ b ≈ D. For example, in the case of ds-DNA in a solution at high salt concentration the distance between subsequent base pairs is 0.34 nm, but D can be chosen as the diameter of the helix, D = 2 nm, so that there are 6 base pairs per subunit. For such conditions, the persistence length is roughly ℓ p = 50 nm so the corresponding stiffness parameter of the model presented in section 2.5 would be κ = ℓ p /ℓ b = 25.
Single semiflexible polymers (as they physically occur in very dilute solutions under good solvent conditions) are typically well described by the Kratky-Porod (KP) worm-like chain model when ℓ p ≫ D. However, apart from extremely long chains, this model is inappropriate for 'thick worms' with D and ℓ p of the same order, which is found sometimes for bottlebrush polymers, dendronized polymers, etc. Also semiflexible polymers strongly adsorbed on attractive substrates have to be considered with care: their effective persistence length, measured from the decay of bond-orientational correlations along the contour of the chain molecule, exhibits crossovers from ℓ p to 2ℓ p for distances s along the contour less than 2ℓ p ; for large distances, the power-law decay of orientational correlations implies that the KP model fails qualitatively due to its neglect of excluded volume interactions.
The excluded volume interaction is also the central mechanism for the onset of nematic order in semidilute or concentrated lyotropic solutions of semiflexible polymers. We have discussed various theoretical predictions for the resulting phase diagrams (in terms of the variables osmotic pressure and density of effective monomeric units), focusing both on the classical Onsager-style Khokhlov-Semenov-Odijk-Chen theory and more recent DFT treatments due to the authors of this review. There, excluded volume between interacting chains is described more realistically, and the strength of this excluded volume is appropriately renormalized. The MD work reviewed here has validated this DFT treatment, apart from deviations in the nematic phase, where DFT predicts too much order, due to the neglect of collective 'deflection fluctuations' . However, important problems about bulk nematic phases that are still incompletely understood concern the Frank elastic constants, and how they show up in suitable experimentally observable response functions. So far, DFT includes orientational distribution functions only, without proper coupling between local density and orientational fluctuations; fully conclusive MD work would require the study of considerably larger systems and strongly enhanced statistical effort, which is still beyond reach at the present time.
We have also mentioned the extension of the description to solvents of variable quality, where unmixing in the isotropic phase competes with the isotropic-nematic transition. Finally, effects of confinement in pores with slit geometry have been discussed. As expected, slit pores favor chains parallel to the walls of the slit, and hence lead to the phenomenon of 'capillary nematization' . Further, attractive monomer-surface interactions can give rise to interesting textures of locally ordered structures in wall-attached domains of semiflexible chains. Clearly, a very rich behavior then is possible, and only selected facets of this behavior have been studied and presented here.
Finally, we emphasize that liquid crystalline polymers offer ample opportunities for the study of defects in liquid crystalline order [146][147][148][149][150], but this aspect has not been covered in the present review for brevity. We have also not considered the response of polymeric liquid crystalline systems to external perturbations, such as flow [199][200][201][202], electromagnetic fields [203], etc. These phenomena are also important for applications, and constitute a promising area for future simulation studies, with the ever increasing computational power of the available computing resources.