Coarsegrained models of complex fluids at equilibrium and unter shear
Coarsegrained models of complex fluids at equilibrium and unter shear
Friederike Schmid
Fakultät für Physik, Universität Bielefeld,
Universitätsstraße 25, D33615 Bielefeld
Email: schmidphysik.unibielefeld.de
Complex fluids exhibit structure on a wide range of length and time scales, and hierarchical approaches are necessary to investigate all facets of their often unusual properties. The study of idealized coarsegrained models at different levels of coarsegraining can provide insight into generic structures and basic dynamical processes at equilibrium and nonequilibrium.
In the first part of this lecture, some popular coarsegrained models for membranes and membrane systems are reviewed. Special focus is given to beadspring models with different solvent representations, and to randominterface models. Selected examples of simulations at the molecular and the mesoscopic level are presented, and it is shown how simulations of molecular coarsegrained models can bridge between different levels.
The second part addresses simulation methods for complex fluids under shear. After a brief introduction into the phenomenology (in particular for liquid crystals), different nonequilibrium molecular dynamics (NEMD) methods are introduced and compared to one another. Application examples include the behavior of liquid crystal interfaces and lamellar surfactant phases under shear. Finally, mesoscopic simulation approaches for liquid crystals under shear are briefly discussed.
1 Introduction
The term “Complex fluid” or “soft condensed matter” – these two are often used interchangeably – refers to a broad class of materials, which are usually made of large organic molecules and have a number of common features [2, 3]: They display structure on a nanoscopic or mesoscopic scale; characteristic energies are of the order of at room temperature, hence the properties are to a large extent dominated by entropic effects, and the materials respond strongly to weak external forces [4]; the characteristic response times span one or several orders of magnitude, and the rheological properties of the fluids are typically nonNewtonian [3, 5]. Some examples are polymer melts and solutions [6, 7, 8, 9], emulsions, colloids [10, 11, 12], amphiphilic systems [13, 14], and liquid crystals [15, 16, 17].
Computer simulations of complex fluids are particularly challenging, due to the hierarchy of length and time scales that contributes to the material properties. With current computer resources, it is practically impossible to describe all aspects of such a material within one single theoretical model. Therefore, one commonly uses different models for different length and time scales. The idea is to eliminate the smallscale degrees of freedom in the largescale models, and to incorporate the details of the smallscale properties in the parameters of a new, simpler, model. This is called coarsegraining.
There are two aspects to coarsegraining. First, systematic coarsegraining procedures must be developed and applied in order to study materials quantitatively on several length and time scales. This is the challenge of multiscale modeling, one major growth area in materials science. It is, however, not the subject of the present lecture. Some recent reviews can be found in Refs. [18, 19]
Second, coarsegrained idealized models are used to study generic properties of materials on a given scale, and physical phenomena which are characteristic for a whole class of materials. Here, microscopic details are disregarded because they are deemed irrelevant for the physical properties under consideration. This approach to materials science has a long tradition, going back to Ising’s famous lattice model for magnetism.
In the present lecture, we shall concentrate on coarsegrained models of the second type. We shall discuss how such models can be constructed, how they can be studied by computer simulations, and how such studies help to improve our understanding of equilibrium and nonequilibrium properties of complex materials. We will not attempt to present an overview – this would require a separate book – but rather focus on selected case studies. In the first part, we will consider coarsegrained model systems for amphiphilic membranes and their use for studies of equilibrium properties. In the second part, we will turn to complex fluids under shear, i.e., nonequilibrium systems, and discuss coarsegrained approaches for membrane systems and liquid crystals.
2 Coarsegrained models for surfactant layers
In this section, we shall discuss coarsegrained simulation models and methods for amphiphilic membranes and membrane stacks. After a brief introduction into the phenomenology and some general remarks on coarsegraining, we will focus on two particular types of coarsegrained models: Particlebased beadspring models and mesoscopic membrane models. We shall present some typical examples, and illustrate their use with applications from the literature.
2.1 Introduction
Amphiphilic membranes are a particular important class of soft material structures, because of their significance for biology [20, 21, 22]. Biomembranes are omnipresent in all living beings, they delimit cells and create compartments, and participate in almost all biological functions. Figure 1 shows an artist’s view on such a biomembrane. It is mainly made of two coupled layers of lipids, which serve as a matrix for embedded proteins and other molecules.
This structure rests on a propensity of lipids to selfassemble into bilayers when exposed to water. The crucial property of lipids which is responsible for this behavior is their amphiphilic character, i.e., the fact that they contain both hydrophilic (water loving) and hydrophobic (water hating) parts. Most lipid molecules have two non polar (hydrophobic) hydrocarbon chain (tails), which are attached to one polar (hydrophilic) head group. Figure 2) shows some typical examples. The tails vary in their length, and in the number and position of double bonds between carbon atoms. Whereas single bonds in a hydrocarbon chain are highly flexible, in the sense that chains can rotate easily around such a bond, double bonds are stiff and may enforce kinks in the chain. Chains with no double bonds are called saturated. The variety of polar head groups is even larger. Overviews can be found in [20].
By assembling into bilayers, the lipids can arrange themselves such that the head groups shield the tails from the water environment. Therefore, lipids in lipidwater mixtures often tend to form lamellar stacks of membranes separated by thin water layers (Fig. 3). Such stacks can typically hold up to 20 % water. In the more dilute regime, membranes may close up to form vesicles or other more disordered structures. Vesicles play an important role in biological systems, since they provide closed compartments which can be used to store or transport substances. Single, isolated membranes are metastable with respect to dissociation for entropic reasons, but they have very long lifetimes.
Depending on the particular choice of lipids and the properties of the surrounding aqueous fluid (pH, salt content etc.), several other structures can also be observed in lipidwater systems – spherical and cylindrical micelles, ordered structures with hexagonal or cubic symmetry etc. [23]. In this lecture, we shall only regard bilayer structures.
Membranes also undergo internal phase transitions. Figure 4 shows schematically some characteristic phases of onecomponent membranes. The most prominent phase transition is the “main” transition, the transition from the liquid state into one of the gel states, where the layer thickness, the lipid mobility, and the conformational order of the chains, jump discontinuously as a function of temperature. In living organisms, most membranes are maintained in the liquid state. Nevertheless, the main transition presumably has some biological relevance, since it is very close to the body temperature for some of the most common lipids (e.g., saturated phospholipids) [20].
In sum, several aspects of membranes must be studied in order to understand their structure and function, in biological as well as in artificial systems:

Selfassembling mechanisms

The inner structure and internal phase transitions, and the resulting membrane properties (permeability, viscosity, stiffness)

The organization of membranes on a mesoscopic scale (vesicles, stacks, etc.)

The interaction of membranes with macromolecules.
These topics relate to the physics and chemistry of amphiphiles on many length and time scales, and several driving forces contribute to the unique properties of membranes: On the scale of atoms and small molecules, one has the forces which are responsible for the selfassembly: The hydrophobic effect, and the interactions between water and the hydrophilic head groups. On the molecular scale, one has the interplay of local chain conformations and chain packing, which is responsible for the main transition and determines the local structure and fluidity of the membrane. On the supramolecular scale, one has the forces which govern the mesoscopic structure and dynamics of membrane systems: The membrane elasticity, the thermal fluctuations, the hydrodynamic interactions with the surrounding fluid, the thermodynamic forces that drive phase separation in mixed membranes etc.
Studying all these aspects within one single computer simulation model is clearly impossible, and different models have been devised to study membranes at different levels of coarsegraining. Among the oldest approaches of this kind are the Doniach [24] and the Pink model [25, 26], generalized Ising models, which have been used to investigate almost all aspects of membrane physics on the supramolecular scale [22, 27]. With increasing computer power, numerous other simulation models have become treatable in the last decades, which range from atomistic models, molecular coarsegrained models, up to mesoscopic models that are based on continuum theories for membranes.
In the following, we shall focus on two important examples: Coarsegrained chain models designed to describe membranes on the molecular scale, and random interface models that describe membranes on the mesoscopic scale. We stress that this our overview is not complete, and there exist other successful approaches, which are not included in our presentation. A relatively recent review can be found in Ref. [28].
2.2 Coarsegrained molecular models
In coarsegrained molecular models, molecules are still treated as individual particles, but their structure is highly simplified. Only the ingredients deemed essential for the material properties are kept. For amphiphiles, these are the amphiphilic character (for the selfassembly), and the molecular flexibility (for the internal phase transitions).
Springbead models
The first coarse grained molecular models have been formulated on a lattice, relatively shortly after the introduction of the abovementioned Pink model [25, 26]. One particularly popular model is the Larson model [29, 30]: Water molecules () occupy single sites of a cubic lattice, and amphiphile molecules are represented by chains of “tail” () and “head” () monomers. Only particles on neighbor lattice sites interact with each another. The lattice is entirely filled with , , or particles. If one further takes all hydrophilic particles, and , to be identical, the interaction energy is determined by a single interaction parameter, which describes the relative repulsion between hydrophobic and hydrophilic particles. The model reproduces selfassembly and exhibits many experimentally observed phases, i.e., the lamellar phase, the hexagonal phase, the cubic micellar phase, and even the bicontinuous gyroid phase [30, 31].
Lattice models can be simulated efficiently by Monte Carlo methods. However, they have the obvious drawback of imposing an a priori anisotropy on space, which restricts their versatility. For example, internal membrane phase transitions and tilted phases cannot be studied easily. In dynamical studies, one is restricted to using Monte Carlo dynamics, which is unrealistic. Therefore, most more recent investigations rely on offlattice models.
Offlattice molecular approaches to study selfassembling amphiphilic systems have been applied for roughly 15 years. Presumably the first model was introduced by B. Smit et al. [32] in 1990. A schematic sketch is shown in Fig. 5 a). As in the Larson model, the amphiphiles are represented by chains made of very simple  or units, which are in this case spherical beads. “Water” molecules are represented by free beads. Beads in a chain are joined together by harmonic springs, where a cutoff is sometimes imposed in order to ensure that the beads cannot move arbitrarily far apart from one another. Nonbonded beads interact via simple shortrange pairwise potentials, e.g., a truncated LennardJones potential. The parameters of the potentials are chosen such that pairs and pairs effectively repel each other. The Smit model reproduces selfassembly, micelle formation and membrane formation [33, 34, 35, 36, 37].
Most molecular offlattice models are modifications of the Smit model [28, 38, 39, 40]. They differ from one another in the specific choice of the interaction potentials between nonbonded beads, in the bond potentials between adjacent beads on the chain, in the chain architecture (number of head and tail beads, number of tails), and sometimes contain additional potentials such as bending potentials. Nowadays, Smit models are often studied with dissipative particle (DPD) dynamics, and the interaction potentials are the typical soft DPD potentials [38, 40]. A “water particle” is then taken to represent a lump of several water molecules. An introduction into the DPD method is given in B. Dünweg’s lecture (this book). It works very well for membrane simulations as long as the time step is not too large [41].
This approach is presented in more detail in B. Smit’s seminar (this book).
Implicit solvent and phantom solvent
The primary requirement for a lipid membrane model is to ensure that the lipids maintain a bilayer structure. In the models discussed in the previous section 2.2.1, the bilayers are stabilized by the surrounding solvent particles. If the interactions between , , and particles are chosen suitably, the solvent particles and the amphiphiles arrange themselves such that the hydrophilic and hydrophobic particles are separated (microphase separation).
This “explicit solvent” approach is in some sense the most natural one and has several advantages: The bilayers are fully flexible, and they are surrounded by a fluid. The fact that the fluid is composed of large beads as opposed to small water molecules is slightly problematic, since the solvent structure may introduce correlations, e.g., in simulations with periodic boundary conditions. However, the correlations disappear if the membranes or their periodic images are separated by a large amount of solvent.
On the other hand, explicit solvent models also have a disadvantage: The simulation time depends linearly on the total number of beads. A large amount of computing time is therefore spent on the uninteresting solvent particles. Therefore, effort has been spent on designing membrane models without explicit solvent particles.
One early idea has been to force the amphiphiles into sheets by tethering the head groups to twodimensional opposing surfaces, e.g., by a harmonic potential, or by rigid constraints [42, 43, 44, 45, 46] (Fig. 5 b)). This model, which we call “sandwich model”, is simple, extremely efficient, and the configurations are easy to analyze, since the bilayer is always welldefined. On the other hand, the bilayer has no flexibility, i.e., undulations and protrusions [47] are supressed and important physics is lost.
A second approach is to eliminate the solvent, and represent its effect on the amphiphiles by appropriate effective interactions between monomers (Fig. 5 c)). This way of modeling solvent is common in polymer simulations. Nevertheless, producing something as complex as membrane selfassembly is a nontrivial task. The first implicit solvent model was proposed only rather recently, in 2001, by Noguchi and Takasu [48]. They mimick the effect of the solvent by a multibody “hydrophobic potential”, which depends on the local density of hydrophobic particles. From a physical point of view, the introduction of multibody potentials in an implicit solvent model is reasonable – multibody potentials emerge automatically whenever degrees of freedom are integrated out systematically. Nevertheless, people usually tend to favor pair potentials. Farago [49] and Cooke et al. [50] have shown that it is possible to stabilize fluid bilayers with purely pairwise interactions in a solventfree model. The crucial requirement seems to be that the interactions between hydrophobic beads have a smooth attractive part of relatively long range [50].
The clear benefit of implicit solvent models is their high efficiency. They are very well suited to study large scale problems like vesicle formation [48, 50], vesicle fusion [51], vesicles subject to external forces [52], phase separation on vesicles [50] etc. On the other hand, they also have limitations: They require tricky potentials; the pressure is essentially zero (since the membranes are surrounded by empty space); also, they cannot be used for dynamical studies where the hydrodynamical coupling with the solvent becomes important.
A third approach was recently put forward by ourselves: the “phantom solvent” model [46] (Fig. 5d)). We introduce explicit solvent particles, which do not interact with one another, only with the amphiphiles. The amphiphiles perceive them as repulsive soft beads. In the bulk, they simply form an ideal gas. They have the simple physical interpretation that they probe the free volume which is still accessible to the solvent, once the membrane has selfassembled. Thus the selfassembly is to a large extent driven entropically.
In Monte Carlo simulations, this model is just as efficient as solventfree models, because Monte Carlo moves of phantom beads that are not in contact with a membrane are practically free of (computational) charge. It is robust, we do not need to tune the lipid potentials in order to achieve selfassembly. Pressure can be applied if necessary, and with a suitable dynamical model for the solvent (e.g., dissipative particle dynamics), we can also study (hydro)dynamical phenomena. Compared with Smittype models, the phantom solvent model has the advantage that the solvent is structureless and cannot induce artificial correlations. On the other hand, the fact that it is compressible like a gas, may cause problems in certain dynamical studies.
To summarize, we have presented several possible ways to model the solvent environment of a selfassembled membrane, and discussed the advantages and disadvantages. The different models are sketched schematically in Fig. 5.
We shall now illustrate the potential of molecular coarsegrained models with two examples.
First application example: The ripple phase
Our first example shows how coarsegrained simulations can help to understand the inner structure of membranes. We have introduced earlier the liquid and gel phases in membranes (Fig. 4). However, one rather mysterious phase was missing in our list: The ripple phase (). It emerges as an intermediate state between the tilted gel phase () and the liquid phase (), and according to freezeetch micrographs, AFM pictures, and Xray diffraction studies, it is characterized by periodic wavy membrane undulations. More precisely, one observes two different rippled structures: Figure 6 shows electron density maps of these two states, as calculated from Xray data by Sengupta et al [53]. The more commonly reported phase is the asymmetric ripple phase, which is characterized by sawtooth profile and has a wavelength of about 150 Angstrom. The symmetric rippled phase tends to form upon cooling from the liquid phase and has a period twice as long. It is believed to be metastable. The molecular structure of both rippled states has remained unclear until very recently.
Atomistic simulations of de Vries et al, have finally shed light on the microscopic structure of the asymmetric ripple phase [54]. These authors observed an asymmetric ripple in a Lecithin bilayer, which had a structure that was very different from the numerous structures proposed earlier. Most strikingly, the rippled membrane is not a continuous bilayer. The ripple is a membrane defect where a single layer spans the membrane, connecting the two opposing sides (Fig. 7).
The simulations of de Vries et al seem to clarify the structure of the asymmetric ripple phase, at least for the case of Lecithin. Still the question remains whether their result can be generalized to other lipid layers, or whether it is related to the specific properties of Lecithin head groups. This question can be addressed with idealized coarsegrained models.
Rippled phases had also been observed in simulations of coarsegrained models. Kranenburg et al have reported the existence of a rippled state in a Smit model with soft DPD interactions [55, 56]. The structure of their ripple is different from that of Fig. 7. Judging from this result, one is tempted to conclude that the structure of de Vries et al may not be generic.
However, we recently found ripples with a structure very similar to that of Fig. 7 in simulations of a simple beadspring model by (Lenz and Schmid [57, 58]). Our lipids are represented by linear chains with six tail beads (size ) and one head bead (size ), which are connected by springs of size . Tail beads attract one another with a truncated and shifted LennardJones potential. Head beads are slightly larger and purely repulsive. The solvent environment is modelled by a fluid of phantom beads of size . The system was studied with Monte Carlo simulations at constant pressure, using a simulation box of fluctuating size and shape in order to ensure that the pressure tensor is diagonal.
Membranes were found to selfassemble spontaneously at sufficiently low temperatures, they exhibit a fluid phase as well as a tilted gel phase. The two phases are separated by a temperature region where ripples form spontaneously. In small systems, the ripples are always asymmetric. In larger systems, asymmetric ripples form upon heating from the phase. Upon cooling from the phase, a second, symmetric, type of ripple forms, which has twice the period than the asymmetric ripple and some structural similarity, except that the bilayer remains continuous (see Fig. 8). The sideview of this second ripple has striking similarity with the density maps of Fig. 6 b).
The coarsegrained simulations of Lenz thus not only demonstrate that the asymmetric ripple structure found by de Vries is generic, and can be observed in highly idealized membrane models as well. They also suggest a structural model for the yet unresolved symmetric ripple.
Second application example: Membrane stacks
In our second example, we consider the fluctuations and defects in stacks of membranes. We shall show how simulations of coarsegrained membrane models can be used to test and verify mesoscopic models, and how they can bridge between coarsegraining levels.
At the mesoscopic level, membranes are often represented by random interfaces (see also Section 2.3). One of the simplest theories for thermal fluctuations in membrane stacks is the “discrete harmonic model” [59]. It describes membranes without surface tension, and assumes that the fluctuations are governed by two factors: The bending stiffness of single membranes, and a penalty for compressing or swelling the stack, which is characterized by a compressibility modulus . The free energy is given in harmonic approximation
(1) 
where is the average distance between layers. The first part describes the elasticity of single membranes, and the second part accounts for the interactions between membranes. Being quadratic, the free energy functional (1) is simple enough that statistical averages can be calculated exactly.
We have tested this theory in detail with coarsegrained molecular simulations (C. Loison et al [60]) of a binary mixture of amphiphiles and solvent, within a Smittype model originally introduced by Soddemann et al [61]. The elementary units are spheres with a hard core radius , which may have two types: “hydrophilic” or “hydrophobic”. Beads of the same type attract each other at distances less than . “Amphiphiles” are tetramers made of two hydrophilic and two hydrophobic beads, and “solvent” particles are single hydrophilic beads. We studied systems of up to 153600 beads with constant pressure molecular dynamics, using a simulation box of fluctuating size and shape, and a Langevin thermostat to maintain constant temperature.
At suitable pressures and temperatures, the system assumes a fluid lamellar phase. Our configurations contained up to 15 stacked bilayers, which contained 20 volume percent solvent. A configuration snapshot is shown in Fig. 9 (left). In order to test the theory (1), one must first determine the local positions of the membranes in the stack. This was done as follows:

The simulation box was divided into cells. (). The size of the cells fluctuates because the dimensions of the box fluctuate.

In each cell, the relative density of tail beads was calculated. It is defined as the ratio of the number of tail beads and the total number of beads.

The hydrophobic space was defined as the space where the relative density of tail beads is higher than a given threshold . (The value of the threshold was roughly 0.7, i.e., 80 % of the maximum value of ).

Cells belonging to the hydrophobic space are connected to clusters. Two hydrophobic cells that share at least one vortex are attributed to the same cluster. Each cluster defines a membrane.

For each membrane and each position , the two heights and , where the density passes through the threshold , are determined. The mean position is defined as the average
The algorithm identifies membranes even if they have pores. At the presence of other defects, such as necks or passages (connections between membranes), additional steps must be taken, which are not described here. A typical membrane conformation is shown in Fig. 9 (right).
The statistical distribution of can be analyzed and compared with theoretical predictions. For example, the transmembrane structure factor, which describes correlations between membrane positions in different membranes is given by [60]
(2) 
where is the Fourier transform of in the plane and is a dimensionless parameter. The function can also be calculated explicitly [60]. The prediction (2) can be tested by simply plotting the ratio vs. for different . The functional form of the curves should be given by the expression in the square brackets, with only one fit parameter . Figure 10 shows the simulation data. The agreement with the theory is very good over the whole range of wave vectors .
Thus the molecular simulations confirm the validity of the mesoscopic model (1). Moreover, the analysis described above yields a value for the phenomenological parameter . By analyzing other quantities, one can also determine and separately. This gives the elastic parameters of membranes and their effective interactions, and establishes a link between the molecular and the mesoscopic description.
Many important characteristics of membranes not only depend on their elasticity, but also on their defects. For example, the permeability is crucially influenced by the number and properties of membrane pores. A number of atomistic and coarse grained simulation studies have therefore addressed pore formation [62, 63, 64, 65, 66], mostly in membranes under tension. In contrast, the membranes in the simulations of Loison et al have almost zero surface tension. This turns out to affect the characteristics of the pores quite dramatically [67, 68].
Figure 11 shows a snapshot of a hydrophobic layer which contains a number of pores. The pores have nucleated spontaneously. They “live” for a while, grow and shrink without diffusing too much, until they finally disappear. Most pores close very quickly, but a few large ones stay open for a long time. A closer analysis shows that pores are hydrophilic, i.e., the amphiphiles rearrange themselves at the pore edge, such that solvent beads in the pore center are mainly exposed to head beads. The total number of pores is rather high in our system. This is because the amphiphiles are short and the membranes are thin.
The analysis algorithm presented above not only localizes membranes, but also identifies pores, their position and their shape. Therefore one can again test the appropriate mesoscopic theories for pore formation. The simplest Ansatz for the free energy of a pore with the area and the contour length has the form [69]
(3) 
where is a core energy, a material parameter called line tension, and the surface tension. The second term describes the energy penalty at the pore rim. The last term accounts for the reduction of energy due to the release of surface tension in a stretched membrane. In our case, the surface tension is close to zero (), because the simulation was conducted such that the pressure tensor is isotropic. Thus the last term vanishes.
In this simple free energy model, the pore shapes should be distributed according to a Boltzmann distribution, . Figure 12 (left) shows a histogram of contour lengths . The bare data do not reflect the expected exponential behavior. Something is missing. Indeed, a closer look reveals that the naive exponential Ansatz disregards an important effect: The “free energy” (3) gives only local free energy contributions, i.e., those stemming from local interactions and local amphiphile rearrangements. In addition, one must also account for the global entropy of possible contour length conformations. Therefore, we have to evaluate the “degeneracy” of contour lengths , and test the relation
(4) 
Figure 12 (right) demonstrates that this second Ansatz describes the data very well. From the linear fit to the data, one can extract a value for the line tension .
The model (3) makes a second important prediction: Since the free energy only depends on the contour length, pores with the same contour length are equivalent, and the shapes of these pores should be distributed like those of twodimensional selfavoiding ring polymers. In other words, they are not round, but have a fractal structure. From polymer theory, one knows that the size of a two dimensional selfavoiding polymer scales roughly like with the polymer length . In our case, the “polymer length” is the contour length . Thus the area of a pore should scale like
(5) 
Figure 13 shows that this is indeed the case.
Other properties of pores have been investigated, e.g., dynamical properties, pore lifetimes, pore correlations etc. [67]. All results were in good agreement with the line tension model (3). In sum, we find that the fluctuations and (pore) defects in membrane stacks can be described very well by a combination of two simple mesoscopic theories: An effective interface model for membrane undulations (Eq. (1)), and a line tension model for the pores (Eq. (3) [68].
The simulations of Loison et al thus demonstrate nicely that simulations of coarsegrained molecular models can be used to test mesoscopic theories, and to bridge between the microscopic and the mesoscopic level.
2.3 Mesoscopic membrane models
When it comes to large scale structures and long time dynamics, simulations of molecular models become too cumbersome. It then proves useful to drop the notion of particles altogether and work directly with continuous mesoscopic theories. In general, continuum theories are mostly formulated in terms of continuous fields, which stand for local, coarsegrained averages of some appropriate microscopic quantities. In the case of membranes systems, it is often more adequate to keep the membranes as separate objects, despite their microscopic thickness, and treat them as interfaces in space. This leads to randominterface theories.
Random interface models
In randominterface models, membranes are represented by twodimensional sheets in space, whose conformations are distributed according to an effective interface free energy functional. An example of such a functional has already been given in Eq. (1). However, the free energy (1) can only be used to describe slightly fluctuating planar membranes in the plane. In order to treat membranes of arbitrary shape (vesicles etc.), one must resort to the formalism of differential geometry [13].
One of the simplest approaches of this kind is the spontaneous curvature model. The elastic free energy of a membrane is given by [71]
(6) 
where integrates over the surface of the membrane, and and are the mean and Gaussian curvature, which are related to the local curvature radii and via , [70]. The parameters , (the bending stiffnesses) and (the spontaneous curvature) depend on the specific material properties of the membrane. The total membrane surface is taken to be constant.
In the case and , the free energy (6) describes membranes that want to be flat, and imposes a bending penalty on all local curvatures. implies that the membrane has a favorite mean curvature. It is worth noting that the integral over the Gaussian curvature, , depends solely on the topology of the interfaces. For membranes without rims, the GaussBonnet theorem states that
(7) 
where the Euler characteristics counts the number of closed surfaces (including cavities) minus the number of handles .
Other free energy expressions have been proposed, which allow for an asymmetric distribution of lipids between both sides of the membrane, and/or for fluctuations of [72, 73]. Here, we shall only discuss the spontaneous curvature model. Our membranes shall be selfavoiding, i.e., they have excluded volume interactions and cannot cross one another.
Our task is to formulate a simulation model that mimicks the continuous space theory as closely as possible. More precisely, we first need a method that samples random selfavoiding surfaces, and second a way to discretize the free energy functional (6) on those surfaces.
A widely used method to generate selfavoiding surfaces is the dynamic triangulation algorithm [74, 75, 76, 77, 78, 79]. Surfaces are modeled by triangular networks of spherical beads, which are linked by tethers of some given maximum length (tetheredbead model). The tethers define the neighbor relations on the surface. In order to enforce selfavoidance, the beads are equipped with hard cores (diameter ), and the maximum tether length is chosen in the range . The surfaces are sampled with Monte Carlo, with two basic types of moves: Moves that change the position of beads, and moves that change the connectivity of the network, i.e., redistribute the tethers between the chains. The latter is done by randomly cutting tethers and flipping them between two possible diagonals of neighboring triangles (Fig. 14). Updates of bead positions are subject to the constraint that the maximum tether length must not exceed . In addition, connectivity updates have to comply with the requirement that at least three tethers are attached to one bead, and that two beads cannot be connected to one another by more than one tether. Otherwise, attempted updates are accepted or rejected according to the standard Metropolis prescription. Of course, additional, more sophisticated Monte Carlo moves can be introduced as well.
Next, the free energy functional (6) must be adapted for the triangulated surface. If no rims are present, we only need an expression for the first term. Since the integral over the Gaussian curvature is then essentially given by the Eulercharacteristic, it is much more convenient and accurate to evaluate the latter directly (i.e., count closed surfaces minus handles) than to attempt a discretization of .
The question how to discretize the bending free energy in an optimal way is highly nontrivial, and several approximations have been proposed over the years [79]. Here we give a relatively recent expression due to Kumar et al in 2001 [80]: The surface integral is replaced by a sum over all triangles with area . For each triangle , one considers the three neighbor triangles and determines the length of the common side as well as the angle between the two triangle normal vectors. The free energy is then approximated by
(8) 
This completes the tetheredbead representation for random interfaces. We have presented one of the most simple variants. The models can be extended in various directions, e.g., to represent mixed membranes [81], membranes that undergo liquid/gel transitions [82], membranes with fluctuating geometry [83] or pores [84] etc. They have even been used to study vesicle dynamics, e.g., vesicles in shear flow [85, 86]. The dynamics can be made more realistic by using hybrid algorithms, where the beads are moved according to molecular dynamics or Brownian dynamics, and the tethers are flipped with Monte Carlo. Unfortunately, the tether flips cannot easily be integrated in a pure molecular dynamics simulation. This is a slight drawback of the tetheredbead approach.
An interesting alternative has recently been proposed by Noguchi et al [87]. He developed a specific type of potential, which forces (single) beads to selfassemble into a twodimensional membrane without being connected by tethers. The membrane is coarsegrained as a curved surface as in tetheredbead models, but tethers are no longer necessary. This opens the way for a fully consistent treatment of dynamics in membrane systems with frequent topology changes.
Application example: passage of vesicles through pores
As an example for an application of a tetheredbead model, we shall discuss a problem from the physics of vesicles. Linke, Lipowsky and Gruhn have recently investigated the question whether a vesicle can be forced to cross a narrow pore by osmotic pressure gradients [88, 89]. Vesicle passage through pores plays a key role in many strategies for drug delivery. On passing through the pore, the vesicle undergoes a conformational change, which is expensive and creates a barrier. On the other hand, the vesicle can reduce its energy, if it crosses into a region with lower osmotic pressure. Linke et al studied the question, whether the barrier height can be reduced by a sufficient amount that the pore is crossed spontaneously.
To this end, they considered a vesicle filled with osmotically active particles. The concentration of these particles outside the vesicle is given by on the two sides of the pore. For a vesicle in the process of crossing the pore from side 1 to side 2, the vesicle area on both sides of the membrane is denoted , and the volumes are . Hence the osmotic contribution to the free energy is
(9) 
where is some reference volume. In the Monte Carlo simulations, a series of different crossing stages was sampled (see Fig. 15). This was done by introducing an appropriate reweighting factor, which kept the area within a certain, predetermined range. Several simulations were done with different windows for . The procedure has similarity with umbrella sampling, except that the windows did not overlap. Nevertheless, the extrapolation of the histograms allowed to evaluate the free energy quite accurately, as a function of .
The resulting free energy landscape is shown in Fig. 16, for different osmolarities at the side where the vesicle starts. As expected, one finds a free energy barrier for low . On increasing , the height of the barrier decreases, until it finally disappears. Hence the vesicle can be forced to cross the pore spontaneously, if is chosen sufficiently high. In contrast, changing the bending rigidity has hardly any effect on the height of the potential barrier.
This application shows how mesoscopic models of membrane systems can be used to study the physical basis of a potentially technologically relevant process. A comparable study with a particlebased model would have been much more expensive, and both the setup and the analysis would have been more difficult: One has no welldefined interface, one has only indirect information on the bending stiffness etc. The simplifications in the model clearly help not only to get results more quickly, but also to understand them more systematically.
2.4 Summary
In this section, we have introduced and discussed two important classes of coarsegrained models for membrane systems: Beadspring models for simulations on the coarsegrained molecular level, and random interface models for the coarsegrained mesoscopic level. We have shown that such studies can give insight into generic properties of membranes and into basic physical mechanisms on different scales. We have also shown that it is possible to bridge between different levels, and that simulations of particlebased models can be used to test mesoscopic theories.
3 Liquid crystals and surfactant layers under shear
So far, we have talked about equilibrium simulations. The special properties of complex fluids lead to particularly peculiar behavior at nonequilibrium. In this section, we will discuss coarsegrained simulations of complex fluids under shear. As in the previous section, we will move from the coarsegrained molecular level to the mesoscopic level. However, we will now focus on the simulation methods rather than on the construction of simulation models. As example systems, we will consider surfactant layers and liquid crystals.
3.1 Introduction
Strain rate and shear stress
We begin with recalling some basic fluid mechanics [90]. In the continuum description, the fluid is thought to be divided into fluid elements, which are big enough that microscopic details are washed out, but small enough that each element can still be considered to be homogoneous. On a macroscopic scale, the fluid is described by spatially varying fields, e.g., the velocity or flow field . The velocity gradient in the direction perpendicular to the flow is called the local strain rate.
One of the central quantities in fluid mechanics is the stress tensor . It describes the forces that the surrounding fluid exerts on the surfaces of a fluid element: For a surface with orientation (), the force per area is . The diagonal components of give the longitudinal stress and correspond to normal forces. The offdiagonal components give the shear stress and correspond to tangential forces.
A fluid element with the velocity is thus accelerated according to
(10) 
where is the local mass density. The stress tensor is often divided in three parts.
(11) 
The first term describes the effect of the local pressure. It is always present, even in an equilibrium fluid at rest. The second term, , accounts for elastic restoring forces, if applicable. The last term, , gives the contribution of viscous forces and is commonly called the viscous stress tensor. It vanishes in the absence of velocity gradients, i.e., if . The concrete relation between the stress tensor and the local fields is called constitutive equation and characterizes the specific fluid material.
In standard hydrodynamics, one considers socalled Newtonian fluids, which have no elasticity () and a linear, instantaneous relation between the viscous stress tensor and the velocity gradient,
(12) 
In the case of an incompressible, isotropic, Newtonian fluid, one has only one independent material parameter, the kinematic viscosity , and the constitutive equation reads
(13) 
where the value of the local pressure is determined by the requirement that the incompressibility condition
(14) 
is always fulfilled. Equations (10) in combination with (13) are the famous Navierstokes equations (without external forces).
Complex fluids have internal degrees of freedom, and large characteristic length and time scales. Therefore, they become nonNewtonian already at low shear stress, and standard hydrodynamics does not apply. The relation between the stress tensor and the velocity gradient is commonly not linear, sometimes one has no unique relation at all. If the characteristic time scales are large, one encounters memory effects, e.g., the viscosity may decrease with time under shear stress [3, 5, 90].
These phenomena are often considered in planar Couette geometry: The flow points in the direction, a velocity gradient in the direction is imposed, and the relation between the strain rate and the shear stress is investigated. For Newtonian fluids, the two are proportional. In complex fluids, the differential viscosity may decrease with increasing strain rate (shear thinning), or increase (shear thickening). In some materials, the stressstrain flow curve may even be nonmonotonic. Since such flow curves are mechanically unstable in macroscopic systems, this leads to phase separation, and the system becomes inhomogeneous (Fig. 17) [91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 104, 105]. The separation of a fluid into two states with distinctly different strain rates is called shear banding. As in equilibrium phase transitions, the signature of the coexistence region is a plateau region in the resulting stressstrain flow curve of the inhomogeneous systems. However, one has two important differences to the equilibrium case: First, the location of the plateau cannot be determined by a Maxwell construction, but has to be calculated by solving the full dynamical equations. Second, the plateau is not necessarily horizontal. Even though mechanical stability requires that two coexisting phases must have the same shear stress, the “plateau” may have a slope or even be curved because different coexisting phases may be selected for different shear rates [101, 106] (see Fig. 17).
Nematic liquid crystals
Next we recapitulate some basic facts on liquid crystals. The building blocks of these materials are anisotropic particles, e.g., elongated molecules, associated structures such as wormlike micelles, or even rodlike virusses. A liquid crystal phase is an intermediate state between isotropic liquid and crystalline solid, which is anisotropic (particles are orientationally ordered), but in some sense still fluid (particles are spatially disordered in at least one direction). Some such structures are sketched in Fig. 18. The phase transitions are in some cases triggered by the temperature (thermotropic liquid crystals), and in other cases by the concentration of the anisotropic particles (lyotropic liquid crystals). Lamellar membrane stacks (see previous section) are examples of smectic liquid crystalline structures. In this section, we shall mainly be concerned with nematics, where the particles are aligned along one common direction, but otherwise fluid.
For symmetry reasons, the transition between the isotropic fluid and the nematic fluid must be first order. The nematic order is commonly characterized by the symmetric traceless order tensor
(15) 
where is a unit vector characterizing the orientation of a particle. In an isotropic fluid, is zero. Otherwise, the largest eigenvalue gives the nematic order parameter , and the corresponding eigenvector is the direction of preferred alignment, the “director” . Note that only the direction of matters, not the orientation, i.e., the vectors and are equivalent.
Since the ordered nematic state is degenerate with respect to a continuous symmetry (the direction of the director), it exhibits elasticity: The free energy costs of spatial director variations must vanish if the wavelength of the modulations tends to infinity. At finite wavelength, the system responds to such distortions with elastic restoring forces. Symmetry arguments show that these depend only on three independent material constants [17], the Frank constants , and . The elastic free energy is given by [15]
(16) 
Figure 19 illustrates the corresponding fundamental distortions, the splay (), twist (), and bend mode ().
How does shear affect a liquid crystalline fluid? In the isotropic phase and at low strain rate , the situation is still comparably simple. The velocity gradient imposes a background angular velocity
(17) 
on the fluid, which the particles pick up [107]. Since they rotate faster while perpendicular to the flow, they also acquire a very slight effective orientation at the angle relative to the flow direction. At higher strain rate, the order parameter increases, and the alignment angle decreases [108]. This sometimes even leads to a shearinduced transition into the nematic phase [109, 110, 111, 112, 113, 114, 115].
Deep in the nematic phase, the nematic fluid can be described by the flow field and the director field (). We shall now assume incompressibility () and a linear, instantaneous relation stressstrain relation as in Eq. (12). Compared to simple fluids, there are several complications.
First, nematic fluids are elastic, hence the stress tensor (11) has an elastic contribution
(18) 
The free energy is given by Eq. (16).
Second, the fluid is locally anisotropic, therefore the constitutive equation for the viscous stress tensor contains more independent terms than in simple fluids. The structure of which is compatible with all symmetry requirements reads [16]
(19)  
where is the symmetric flow deformation tensor, , and the vector gives the rate of change of the director relative to the angular velocity of the fluid, . The parameters are the socalled Leslie coefficients. They are linked by one relation [16],
(20) 
hence one has five independent material constants.
Third, the flow field and the director field are coupled. Thus one must take into account the dynamical evolution of the director field,
(21) 
which is driven by the intrinsic body force , and the director stress tensor . Here is the moment of inertia per unit volume. As for the stress tensor, one can derive the general form of the constitutive equations for and with symmetry arguments. The individual expressions can be found in the literature [16]. Here, we only give the final total dynamical equation for ,
(22) 
where , , and is a Lagrange parameter which keeps constant.
To summarize, the dynamical equations of nematohydrodynamics for incompressible nematic fluids are given by Eqs. (10) and (22) together with (11), (16), (18), and (19). These are the LeslieEricksen equations of nematohydrodynamics [16]. The theory depends on five independent viscosity coefficients and on three independent elastic constants . Not surprisingly, it predicts a rich spectrum of phenomena. For example, in the absence of director gradients, steady shear flow is found to have two opposing effects on the director: It rotates the molecules, and it aligns them. The relative strength of these two tendencies depends on the material parameters. As a result, one may either encounter stable flow alignment or a rotating “tumbling” state [116, 117, 118, 119].
The EricksenLeslie theory describes incompressible and fully ordered nematic fluids. The basic version presented here does not account for the possibility of disclination lines and other topological defects, and it does not include the coupling with a variable density and/or order parameter field. The situation is even more complicated in the vicinity of stable or metastable (hidden) thermodynamic phase transitions. Thermodynamic forces may contribute to a nonmonotonic stressstrain relationship, which eventually lead to mechanical phase separation [99].
3.2 Simulating shear on the particle level: NEMD
In the linear response regime, i.e., in the lowshear limit, many rheological properties of materials can be determined from equilibrium simulations. For example, the viscosity coefficients introduced in the previous section can be related to equilibrium fluctuations with GreenKubo relations [120, 121].
However, the nonNewtonian character of most complex fluids and their resulting unique properties manifest themselves only beyond the linear response regime. To study these, nonequilibrium simulations are necessary. We will now briefly review nonequilibrium molecular dynamics (NEMD) methods for the simulation of molecular model systems under shear. Some of the methods are covered in much more detail in Refs. [122, 123, 124].
In NEMD simulations, one faces two challenges. First, one must mechanically impose shear. Second, most methods that enforce shear constantly pump energy into the system. Hence one must get rid of that heat, i.e., apply an appropriate thermostat. We will address these two issues separately.
The most direct way of imposing shear is to confine the system between two rough walls, and either move one of them (for Couette flow), or apply a pressure gradient (for Poiseuille flow). Varnik and Binder have shown that this approach can be used to measure the shear viscosity in polymer melts [125]. It has the merit of being physical: Strain is enforced physically, and the heat can be removed in a physical way by coupling a thermostat to the walls. On the other hand, the system contains two surfaces, and depending on the material under consideration, one may encounter strong surface effects. This can cause problems.
An alternative straightforward way to generate planar Couette flow is to use moving periodic boundary conditions as illustrated in Fig. 21 (LeesEdwards boundary conditions [122, 126]): In order to enforce shear flow with an average strain rate , one proceeds as follows: One replicates the particles in the and direction like in regular periodic boundary conditions,
(23) 
In the direction, the replicated particles acquire an additional velocity and an offset in the direction,
(24) 
Here are the dimensions of the simulation box, and the unit vector in the direction.
LeesEdwards boundary conditions are perfectly sufficient to enforce planar Couette flow. In some applications, it has nevertheless proven useful to supplement them with fictitious bulk forces that favor a linear flow profile. One particularly popular algorithm of this kind is the SLLOD algorithm [123, 127]. For the imposed flow , the SLLOD equations of motion for particles read
(25)  
(26) 
where is the momentum of particle in a reference frame moving with the local flow velocity , and is the regular force acting on the particle . The equations of motion can be integrated with standard techniques, or with specially devised operatorsplit algorithms [129, 130]. It is also possible to formulate SLLOD variants for arbitrary steady flow [128]. The SLLOD algorithm is very useful for the determination of viscosities in complex fluids [131, 120, 132, 133, 134, 135]. It has the advantage that the system relaxes faster towards a steady state, and that the analysis of the data is simplified [123].
A fourth method to enforce shear has recently been proposed by MüllerPlathe [136]. One uses regular periodic boundary conditions, and divides the system into slabs in the desired direction of flow gradient. In periodic intervals, the particle in the central slab () with the largest velocity component in the direction, and the particle in the topmost slab with the largest velocity component in the direction is determined, and the momentum components of the two particles are swapped. This leads to a zigzagshaped shear profile (Fig. 22).
The MüllerPlathe algorithm can be implemented very easily. Just one small change turns a molecular dynamics program for thermal equilibrium into a NEMD program that produces shear flow. Furthermore, the algorithm conserves momentum and energy: Contrary to the algorithms discussed before, it does not pump heat into the system. Instead, the heat is constantly redistributed. It is drained out of the system at and by the momentum swaps, but it is also produced in the bulk through the energy dissipation associated with the shear flow. As a result, the zigzag shear profile is accompanied with a Wshaped temperature profile. Since this is not always wanted, the algorithm is sometimes used in combination with a thermostat. A slight drawback of the algorithm is the fact that it produces intrinsically inhomogeneous profiles, due to the presence of the two special slabs at and . This may cause problems in small systems.
Having introduced these four methods to impose shear, we now turn to the problem of thermostatting. All shearing methods except for the MüllerPlathe method produce heat, therefore the thermostat is an essential part of the simulation algorithm.
A thermostat is a piece of algorithm that manipulates the particle velocities such that the temperature is adjusted to its desired value. This is commonly done by adjusting the kinetic energy. In a flowing fluid with flow velocity , the kinetic energy has a flow contribution and a thermal contribution. Therefore, it is important to define the thermostat in a frame that moves with the fluid, i.e., to couple it to the “peculiar” velocities , rather than to the absolute velocities . Unfortunately, most thermostats don’t account for this automatically, and one has to put in the flow profile by hand. The only exception is the dissipative particle dynamics (DPD) thermostat (see below).
This raises the question how to determine the flow . In a homogeneous system, at low strain rates, it can be acceptable to use the preimposed profile ( in our earlier examples). Thermostats that use the preimposed flow profile are called “profilebiased thermostats”. At high strain rates, and in inhomogeneous systems, the use of profilebiased thermostats is dangerous and may twist the results [123, 134]. In that case, the profiles must be calculated directly from the simulations, “on the fly”. Such thermostats are called “profileunbiased thermostats”.
One can distinguish between deterministic and stochastic thermostatting methods. The simplest deterministic thermostat simply rescales the peculiar velocities of all particles after every time step such that the total thermal kinetic energy remains constant. Among the more sophisticated deterministic thermostats, the ones that have been used most widely in NEMD simulations are the socalled Gaussian thermostats [123], and the Nosé Hoover thermostats [122, 137].
The idea of Gaussian thermostats is very simple: For a system of particles subject to forces , one minimizes the quantity with respect to , with the constraint that either the time derivative of the total energy, or the time derivative of the total kinetic energy be zero (Gaussian isoenergetic and Gaussian isothermal thermostat, respectively). In the absence of constraints, the minimization procedure would yield Newton’s law. The constraints modify the equations of motion in a way that can be considered in some sense minimal. One obtains
(27) 
where is a Lagrange multiplier, which has to be chosen such that the desired constraint – constant energy or constant kinetic energy – is fulfilled.
The Nosé Hoover thermostats are wellknown from equilibrium molecular dynamics and shall not be discussed in detail here. We only note that in NEMD simulations, one uses not only the standard isothermal NoséHoover thermostats, but also isoenergetic variants.
Deterministic thermostats have the slightly disturbing feature that they change the dynamical equations in a rather unphysical way. It is sometimes hard to say how that might affect the system. At equilibrium, NoséHoover thermostats rest on a wellestablished theoretical basis: They can be derived from extended Hamiltonians, they produce the correct thermal averages for static quantities etc. Outside of equilibrium, the situation is much more vague. For the case of steady states, some results on the equivalence of thermostats are fortunately available. Evans and Sarman have shown that steady state averages and time correlation functions are identical for Gaussian isothermal and isoenergetic thermostats and for NoséHoover thermostats [124, 138]. Ruelle has recently proved that the Gaussian isothermal and the Gaussian isoenergetic thermostat are equivalent in an infinite system which is ergodic under space translations [139].
An alternative approach to thermostatting are the stochastic methods, such as the Langevin thermostat [140] and the DPD thermostat [141]. They maintain the desired temperature by introducing friction and stochastic forces. This approach has the virtue of being physically motivated. The dissipative and stochastic terms stand for microscopic degrees of freedom, which have been integrated out in the coarsegrained model, but can still absorb and carry heat. From the theory of coarsegraining, it is wellknown that integrating out degrees of freedom effectively turns deterministic equations of motion into stochastic equations of motion [142, 143, 144].
The Langevin thermostat is the simplest stochastic dynamical model. One modifies the equation of motion for the peculiar momenta by introducing a dissipative friction and a random force ,
(28) 
where fulfills
(29)  
(30) 
i.e., the noise is deltacorrelated and Gaussian distributed with mean zero and variance (fluctuationdissipation theorem). In practice, this can be done by adding and a random number to each force component in each time step . The random number must be distributed with mean zero and variance . Because of the central limit theorem, the distribution does not necessarily have to be Gaussian, if the time step is sufficiently small [145].
The DPD thermostat is an application of the dissipative particle dynamics algorithm [141]. The idea is similar to that of the Langevin algorithm, but instead of coupling to absolute velocities, the thermostat couples to velocity differences between neighbor particles. Therefore, the algorithm is Galilean invariant, and accounts automatically for the difference between flow velocity and thermal velocity. We refer to B. Dünweg’s lecture (this book) for an introduction into the DPD method (see also [137]).
This completes our brief introduction into NEMD methods for systems under shear. We close with a general comment: The strain rate has the dimension of inverse time. It introduces a time scale in the system, which diverges in the limit . Therefore, systems at low strain rates converge very slowly towards their final steady state, and the study of systems at low strain rates is very timeconsuming. In fact, even in coarsegrained molecular models, simulated strain rates are usually much higher than experimentally accessible strain rates.
We now give two examples of recent largescale NEMD studies of inhomogeneous complex systems under shear.
First application example: NematicIsotropic Interfaces under shear
Our first example is a study of the behavior of a NematicIsotropic interface under shear (Germano and Schmid [146, 147]). Interfaces play a central role for shear banding, and understanding their structure should provide a key to the understanding of phase separation under shear. Our study is a first step in this direction. We will discuss the questions whether the structure of an equilibrium interface is affected by shear, whether the phase transition is shifted, and whether shear banding can be observed.
The model system was a fluid of soft repulsive ellipsoids with the aspect ratio 15:1, in a simulation box of size (150:300:150) (in units of the ellipsoid diameter ). The density was chosen in the coexistence region between the nematic and isotropic phase, such that the system phase separates at equilibrium. The initial configuration was a relaxed equilibrium configuration 2.2with a nematic slab and an isotropic slab, separated by two interfaces. The latter align the director of the nematic phase such that it is parallel to the interface [148, 149, 150].
Shear was imposed with LeesEdwards boundary conditions in combination with a profile unbiased Nosé Hoover thermostat, in the direction normal to the interface. The two coexisting phases thus share the same shear stress. Two setups are possible within this “common stress” geometry. In the flowaligned setup, the director points in the direction of flow (the direction); in the logrolling setup, it points in the direction of “vorticity” (the direction). Only the flowaligned setup turned out to be stable in our system. We considered mainly the strain rate . Here is the natural time unit, , with the particle mass , the particle diameter , and the temperature . This strain rate is small enough that the interface is still stable. A configuration snapshot is shown in Fig. 23.
To analyze the system, it was split into columns of size . The columns were further split into bins, which contained enough particles that a local order parameter could be determined. In this manner, one obtains local order parameter profiles for each column in each configuration. From the profile , we determined the local positions of the two interfaces as follows: We computed at least two coarsegrained profiles by convoluting the profile with a symmetrical boxlike coarsegraining function of width . The coarsegraining width must be chosen such that the coarsegrained profiles still reflect the interfaces, but short ranged fluctuations are averaged out. The intersection of two averaged profiles locates a “dividing surface”, where the negative order parameter excess on the nematic side balances the positive order parameter excess on the paranematic side. We define this to be the “local position” of the interface. Due to fluctuations, it depends slightly on the particular choice of the . Nevertheless, the procedure gives remarkably unambiguous values even for strongly fluctuating profiles [150]. The procedure is illustrated in Fig. 24.
Once the positions of the two interfaces have been determined, we calculate profiles for all quantities of interest and shift them by the amount or , respectively. This allows to perform averages over local profiles. Furthermore, the interface positions and themselves can be used to analyze the fluctuation spectrum of the interface positions.
Unexpectedly, the shear turned out to have almost no influence at all on the interfacial structure. Under shear, the interfaces are slightly broadened, indicating that the interfacial tension might be slightly reduced (Fig. 25). But the effect is barely noticeable. Likewise, the density profile and the capillary wave spectrum are almost undistinguishable from those of the equilibrium interface, within the error. The velocity profiles, however, are much more interesting. The flow profile exhibits a distinct kink at the interface position (Fig. 26). Hence we clearly observe shear banding – the total strain is distributed nonuniformly between the two phases. At the same shear stress, the nematic phase supports a higher strain rate than the isotropic phase. Second and surprisingly, the interface apparently also induces a streaming velocity gradient in the vorticity direction (the direction). As a consequence, vorticity flow is generated in opposite directions in the isotropic and in the nematic phase. This flow is symmetry breaking, and as yet, we have no explanation for it. We hope that future theoretical and simulation studies will clarify this phenomenon.
At high strain rates, the interface gets destroyed. Figure 27 shows a nonequilibrium phase diagram, which has been obtained from interface simulations of small systems. The local strain rate is always higher in the nematic phase than in the isotropic phase. It is remarkable that the coexistence region does not close up. Instead, the interface disappears abruptly at average strain rates above .
Second application example: Shearinduced phenomena in surfactant layers
As a second example, we discuss the effect of shear on lamellar stacks of surfactants. This has been studied in coarsegrained molecular simulations by Soddemann, Guo, Kremer et al [151, 152]. The model system was very similar to that used in Section 2.2.4, except that it does not contain solvent particles [61]. Shear was induced with the MüllerPlathe algorithm [136] in combination with a DPD thermostat [61].
The first series of simulations by Guo et al [151] addressed the question, whether the surfactant layers can be forced to reorient under shear. Three different orientations of the layers with respect to the shear geometry were considered: In the transverse orientation, the layers are normal to the flow direction, in the parallel orientation, they are normal to the direction of the velocity gradient, and in the perpendicular orientation, they lie in the plane of the flow and the velocity gradient. The transverse orientation is unstable under shear flow. Experimentally, both a transition from the transverse to the parallel orientation and from the transverse to the perpendicular orientation have been observed [153, 154].
Guo et al have studied this by simulations of a system of dimer amphiphiles. Figures 28 and 29 show a series of snapshots for different times at two different strain rates. Both systems were set up in the transverse state. At low strain rate, the shear generates defects in the transverse lamellae. Mediated by these defects, the lamellae gradually reorient into the parallel state. At high strain rate, the lamellae are first completely destroyed and then reorganize in the perpendicular state.
The final state is not necessarily the most favorable steady state. In fact, Guo et al found that the strain energy dissipation was always smallest in the perpendicular state, regardless of the strain rate. When starting from the transverse state at low strain rates, the parallel state formed nevertheless for kinetic reasons. At intermediate strain rates, shear may induce undulations in the parallel state. This has been predicted by Auernhammer et al [155, 156] and studied by Soddemann et al by computer simulation of a system of layered tetramers [152]. The phenomenon results from the coupling between the layer normal, the tilt of the molecules, and the shear flow. It is triggered by the fact that shear flow induces tilt. The tilted surfactant layer dilates and uses more area, which eventually leads to an undulation instability. This could indeed be observed in the simulations. Figure 30 shows a snapshot of an undulated configuration, and a plot of the undulation amplitude as a function of strain rate. The undulations set in at a welldefined strainrate.
Our two examples show that coarsegrained molecular simulations of complex, inhomogeneous fluids under shear are now becoming possible. A wealth of new, intriguing physics can be expected from this field in the future.
3.3 Simulations at the Mesoscopic Level
We will only touch very briefly on the wide field of mesoscopic simulations for complex fluids under shear.
The challenge of mesoscopic simulations is to find and/or formulate the appropriate mesoscopic model, and then put it on the computer. For liquid crystals, we have already discussed one candidate, the LeslieEricksen theory. However, this is by far not the only available mesoscopic theory for liquid crystals. A popular alternative for the description of lyotropic liquid crystals is the Doi theory, which is based on a Smoluchowski equation for the distribution function for rods. Numerous variants and extensions of the Doi model have been proposed. A relatively recent review on both the LeslieEricksen and the Doi theory can be found in Ref. [157].
The phenomenological equations can be solved with traditional methods of computational fluid dynamics. An alternative approach which is particularly well suited to simulate flows in complex and fluctuating geometries is the Lattice Boltzmann method, which has been discussed by B. Dünweg in his lecture. Yeomans and coworkers have recently proposed a Lattice Boltzmann algorithm for liquidcrystal hydrodynamics, which allows to simulate liquid crystal flow [158, 159]. The starting point is the BerisEdwards theory [160], a dynamical model for liquid crystals that consists of coupled dynamical equations for the velocity profiles and the order tensor (Eq. (15). The LatticeBoltzmann scheme works with the usual discrete velocities and distributions for the fluid flow field on the lattice site . In addition, a second distribution is introduced, which describes the flow of the order tensor field. The time evolution includes free streaming and collision steps. The exact equations are quite complicated and shall not be given here, they can be found in Refs. [159].
As an application example of this method, we discuss a recent simulation by Marenduzzo et al of a sheared hybrid aligned nematic (HAN) cell [161, 162]. The surfaces in a HAN cell impose conflicting orientations on the director of an adjacent nematic fluid, i.e., one surface orients (“anchors”) the fluid in a parallel way (“planar”), and the other in a perpendicular way (“homeotropic”). In the system under consideration, the parameters were chosen such that the fluid is in the isotropic phase, but very close to the nematic phase. The fluid is at rest at time . At times , the surfaces are moved relative to each other, and drag the fluid with them (noslip boundary conditions). Figure 31 shows snapshots of the cell for different times. The evolution of the structure inside the cell is quite complex. First, two ordered bands form close to the surface, due to the fact that the flow field still builds up and strain rates are quite high in the vicinity of the surfaces. The director in these bands is flowaligned, and in the lower band, its direction competes with the anchoring energy of the homeotropic surface. As a result, the lower band unbinds and crosses the system to join the other band. Finally, the top band also unbinds and moves towards the center of the cell. This simulation shows nicely how a complex dynamical process, which results from an interplay of shear flow and elastic deformations, can be studied with a mesoscopic simulation method. We note that the time scale of this simulation is seconds, i.e., inaccessible for molecular simulations.
3.4 Summary
Under shear, complex fluids exhibit a wealth of new, interesting, and practically relevant phenomena. To study these, various simulation approaches for different levels of coarsegraining have been developed. In this section, we have presented and discussed several variants of nonequilibrium molecular dynamics (NEMD), and illustrated their use with examples of large scale NEMD simulations of inhomogeneous liquid crystals and surfactant layers. Furthermore, we have briefly touched on mesoscopic simulation methods for liquid crystals under shear.
4 Conclusions
Due to the rapidly improving performance of modern computers, complex fluids can be studied on a much higher level today than just ten years ago. A large number of coarsegrained models and methods have been designed, which allow to investigate different aspects of these materials, at equilibrium as well as far from equilibrium. The new computational possibilities have vitally contributed to the current boost of soft matter science.
In the two central sections of this lecture, we have given an introduction into two important chapters of computational soft matter: In Section 2, we have given a rough overview over idealized coarsegrained simulation models for membranes, and hopefully conveyed an idea of the potential of such models. In Section 3, we have discussed simulation methods for complex fluids under shear. Generally, nonequilibrium studies of complex fluids are still much more scarce than equilibrium studies. Much remains to be done, both from the point of view of method development and applications, and many exciting developments can be expected for the near future.
Acknowledgments
I would like to thank Olaf Lenz, Claire Loison, Michel Mareschal, Kurt Kremer, and Guido Germano, for enjoyable collaborations that have lead to some of the results discussed in this paper. I am grateful to Alex de Vries, SiewertJan Marrink, Gunnar Linke, and Thomas Gruhn, for allowing to use their configuration snapshots. This work was supported by the John von Neumann institute for computing, and by the Deutsche Forschungsgemeinschaft.
References
 M. Daoud and C. E. Williams, Eds. (1995) Soft Matter Physics. SpringerVerlag, Berlin
 R. G. Larson (1999) The Structure and Rheology of Complex Fluids. Oxford University Press, New York
 J. N. Israelachvili (1992) Intermolecular and Surface Forces. Academic Press, London
 M. Kröger (2003) Simple models for complex nonequilibrium fluids. Physics reports 390, p. 453
 P. J. Flory (1969) Statistical Mechanics of Chain Molecules. Interscience Publishers, New York
 P.G. de Gennes (1979) Scaling Concepts in Polymer Physics. Cornell University Press, Ithaca
 M. Doi (1992) Introduction to Polymer Physics. Clarendon Press, Oxford
 M. Doi and S. Edwards (1986) The Theory of Polymer Dynamics. Clarendon Press, Oxford
 W. B. Russel, D. A. Saville, and W. R. Schowalter (1989) Colloidal Dispersions. Cambridge University Press, Cambridge
 R. J. Hunter (1989) Foundations of Colloid Science. Clarendon Press, Oxford
 M. Borowko, Ed. (2000) Computational Methods in Surface and Colloid Science. Marcel Dekker Inc., New York
 S. A. Safran (1994) Statistical Thermodynamics of Surfaces, Interfaces, and Membranes. AddisonWesley, Reading, MA
 G. Gompper and M. Schick (1994) SelfAssembling Amphiphilic Systems, in Phase Transitions and Critical Phenomena, C. Domb and J. Lebowitz Eds., Academic Press, London, 1994 16, p. 1.
 P. G. de Gennes and J. Prost (1995) The Physics of Liquid Crystals. Oxford University Press, Oxford
 S. Chandrasekhar (1992) Liquid Crystals. Cambridge University Press, Cambridge
 P. M. Chaikin and T. C. Lubensky (1995) Principles of Condensed Matter Physics. Cambridge University Press, Cambridge
 F. MüllerPlathe (2002) Coarsegraining in polymer simulation: From the atomistic to the mesoscopic scale and back. ChemPhysChem 3, p. 754
 S. O. Nielsen, C. F. Lopez, G. Srinivas, and M. L. Klein (2004) Coarse grain models and the computer simulation of soft materials. J. Phys.: Cond. Matter 16, p. R481
 R. B. Gennis (1989) Biomembranes. Springer Verlag, New York
 R. Lipowsky and E. Sackmann Eds. (1995) Structure and Dynamics of Membranes – from Cells to Vesicles. Vol. 1 of Handbook of Biological Physics. Elsevier, Amsterdam
 O. G. Mouritsen (2005) Life – As a Matter of Fat. Springer, Berlin Heidelberg
 K. V. Schubert and E. W. Kaler (1996) Nonionic microemulsions. Ber. Bunseng. Phys. Chemie 100, p. 190
 S. Doniach (1978) Thermodynamic fluctuations in phospholipid bilayers. J. Chem. Phys. 68, p. 4912
 D. A. Pink, T. J. Green, and D. Chapman (1980) Ramanscattering in bilayers of saturated phosphatidylcholines  experiment and theory. Biochemistry 19, p. 349
 A. Caillé, D. A. Pink, F. de Verteuil, and M. Zuckermann (1980) Theoretical models for quasitwodimensional mesomorphic monolayers and membrane bilayers. Can. J. Physique 58, p. 581
 B. Dammann, H. C. Fogedby, J. H. Ipsen, C. Jeppesen, K. Jorgensen, O. G. Mouritsen, J. Risbo, M. C. Sabra, M. M. Sperooto, and M. J. Zuckermann (1995) in Handbook of nonmedical applications of liposomes Vol. 1, Ed. D. D. Lasic, Y. Barenholz, CRC press, p. 85.
 F. Schmid (2000) Systems involving surfactants. In Computational methods in surface and colloid science, Surfactant Science Series Vol. 89, Ed. M. Borowko, Marcel Dekker Inc., New York, p. 631
 R. G. Larson, L. E. Scriven, and H. T. Davis (1985) MonteCarlo simulation of model amphiphilic oilwater systems. J. Chem. Phys. 83, p. 2411
 T. B. Liverpool (1996) Larson models of Amphiphiles in Complex Fluids. In Ann. Rev. Comp. Phys. IV, Ed. D. Stauffer, World Scientific, Singapore, p. 317.
 R. G. Larson (1996) Monte Carlo simulations of the phase behavior of surfactant solutions. J. Physique II 6, p. 1441
 B. Smit, A. G. Schlijper, L. A. M. Rupert, and N. M. van Os (1990) Effects of chainlength of surfactants on the interfacial tension  moleculardynamics simulations and experiments. J. Phys. Chem. 94, p. 6933
 B. Smit, K. Esselink, P. A. J. Hilbers, N. M. van Os, L. A. M. Rupert, and I. Szleifer (1993) Computer simulations of surfactant selfassembly. Langmuir 9, p. 9
 S. Karaborni, K. Esselink, P. A. J. Hilbers, B. Smit, J. Karthäuser, N. M. van Os, and R. Zana, Simulating the selfassembly of gemini (dimeric) surfactants. Science 266, p. 254
 B. J. Palmer and J. Liu (1996) Simulations of micelle selfassembly in surfactant solutions. Langmuir 12, p. 746
 B. J. Palmer and J. Liu (1996) Effects of solutesurfactant interactions on micelle formation in surfactant solutions Langmuir 12, p. 6015
 R. Goetz and R. Lipowsky (1998) Computer simulations of bilayer membranes: Selfassembly and interfacial tension. J. Chem. Phys. 108, p. 7397
 J. C. Shillcock and R. Lipowsky (2002) Equilibrium structure and lateral stress distribution of amphiphilic bilayers from dissipative particle dynamics simulations. J. Chem. Phys. 117, p. 5048
 M. J. Stevens (2004) Coarsegrained simulations of lipid bilayers . J. Chem. Phys. 121, p. 11942
 M. Kranenburg, J. P. Nicolas, and B. Smit (2004) Comparison of mesoscopic phospholipidwater models. Phys. Chem. Chem. Phys. 6, p. 4142
 A. F. Jakobsen, O. G. Mouritsen, and G. Besold (2005) Artifacts in dynamical simulations of coarsegrained model lipid bilayers. J. Chem. Phys. 122, p. 204901
 D. Harries and A. BenShaul (1997) Conformational chain statistics in a model lipid bilayer: Comparison between mean field and Monte Carlo calculations. J. Chem. Phys. 106, p. 1609
 A. Baumgärtner (1995) Asymmetric partitioning of a polymer into a curved membrane. J. Chem. Phys. 103, p. 10669
 A. Baumgärtner (1996) Insertion and hairpin formation of membrane proteins: A Monte Carlo study. Biophys. J. 71, p. 1248
 T. Sintes and A. Baumgärtner (1997) Shortrange attractions between two colloids in a lipid monolayer. Biophys. J. 73, p. 2251
 O. Lenz and F. Schmid (2004) A simple computer model for liquid lipid bilayers. J. Mol. Liquids 117, p. 147
 R. Goetz, G. Gompper, and R. Lipowsky, Mobility and elasticity of selfassembled membranes. Phys. Rev. Lett. 82, p. 221
 H. Noguchi and M. Takasu (2001) Selfassembly of amphiphiles into vesicles: A Brownian dynamics simulation. Phys. Rev. E 64, p. 041913
 O. Farago (2003) ”Waterfree” computer model for fluid bilayer membranes. J. Chem. Phys. 119, p. 596
 I. R. Cooke, K. Kremer, and M. Deserno (2005) Tunable generic model for fluid bilayer membranes. Phys. Rev. E 72, p. 011506
 H. Noguchi and M. Takasu (2001) Fusion pathways of vesicles: A Brownian dynamics simulation. J. Chem. Phys. 115, p. 9547
 H. Noguchi and M. Takasu (2002) Structural changes of pulled vesicles: A Brownian dynamics simulation. Phys. Rev. E 65, p. 051907
 K. Sengupta, V. A. Raghunathan, and J. Katsaras (2003), Structure of the ripple phase of phospholipid multibilayers. Phys. Rev. E 68, p. 031710
 A. H. de Vries, S. Yefimov, A. E. Mark, and S. J. Marrink (2005) Molecular structure of the lecithin ripple phase. PNAS 102, p. 5392
 M. Kranenburg, C. Laforge, and B. Smit (2004) Mesoscopic simulations of phase transitions in lipid bilayers. Phys. Chem. Chem. Phys. 6, p. 4531
 M. Kranenburg and B. Smit (2005) Phase behavior of model lipid bilayers. J. Phys. Chem. B 109, p. 6553
 O. Lenz and F. Schmid (2007) Structure of symmetric and asymmetric ‘ripple’ phases in lipid bilayers. Phys. Rev. Lett. 98, p. 058104
 F. Schmid, D. Düchs, O. Lenz, and B. West (2007) A generic model for lipid monolayers, bilayers, and membranes. Comp. Phys. Comm. 177, p. 168
 N. Lei, C. R. Safinya, and R. F. Bruinsma (1995) Discrete harmonic model for stacked membranes  theory and experiment. J. Phys. II 5, p. 1155
 C. Loison, M. Mareschal, K. Kremer, and F. Schmid (2003) Thermal fluctuations in a lamellar phase of a binary amphiphilesolvent mixture: A molecular dynamics study. J. Chem. Phys. 119, p. 13138
 T. Soddemann, B. Dünweg, and K. Kremer (2001) A generic computer model for amphiphilic systems. Eur. Phys. J. E 6, p. 409
 M. Müller and M. Schick (1996) New mechanism of membrane fusion. J. Chem. Phys. 116, p. 2342
 S.J. Marrink, F. Jähning, and H. Berendsen (1996) Proton transport across transient singlefile water pores in a lipid membrane studied by molecular dynamics simulations. Biophys. J. 71, p. 632
 D. Zahn and J. Brickmann (2002) Molecular Dynamics Study of Water Pores in a Phospholipid Bilayer. Chem. Phys. Lett. 352, p. 441
 T. V. Tolpekina, W. K. den Otter, and W. J. Briels (2004) Simulations of stable pores in membranes: System size dependence and line tension. J. Chem. Phys. 121, p. 8014
 Z. J. Wang and D. Frenkel (2005) Pore nucleation in mechanically stretched bilayer membranes. J. Chem. Phys. 123, p. 154701
 C. Loison, M. Mareschal, and F. Schmid (2004) Pores in bilayer membranes of amphiphilic molecules: CoarseGrained Molecular Dynamics Simulations Compared with Simple Mesoscopic Models. J. Chem. Phys. 121, p. 1890
 C. Loison, M. Mareschal, and F. Schmid (2005) Fluctuations and defects in lamellar stacks of amphiphilic bilayers. Comp. Phys. Comm. 169, p. 99
 J. D. Lister (1975) Stability of lipid bilayers and red bloodcell membranes. Physics Lett. A A 53, p. 193
 E. W. Weisstein (2003) CRC Concise encyclopaedia of mathematics. Chapman & Hall CRC, http://mathworld.wolfram.com
 W. Helfrich (1973) Elastic properties of lipid bilayers  theory and possible experiments. Z. Naturforschung C 28, p. 693
 E. Evans (1974) Bending resistance and chemically induced moments in membrane bilayers. Biophys. J. 14, p. 923
 U. Seifert (1997) Configurations of fluid membranes and vesicles. Adv. Phys. 46, p. 13
 V. A. Kazakov, I. K. Kostov, and A. A. Migdal (1985) Critical properties of randomly triangulated planar random surfaces. Phys. Lett. B 157, p. 295
 A. Billoire and F. David (1986) Scaling properties of randomly triangulated planar random surfaces  a numerical study. Nucl. Phys. B 275, p. 617
 Y. Kantor, M. Kardar, and D. R. Nelson (1986) Statistical mechanics of tethered surfaces. Phys. Rev. Lett. 57, p. 791
 J.S Ho, A. Baumgärtner (1990) Simulations of fluid selfavoiding membranes. Europhys. Lett. 12, p. 295
 D. M. Kroll and G. Gompper (1992) The conformations of fluid membranes  MonteCarlo simulations. Science 255, p. 968
 G. Gompper and D. M. Kroll (1997) Network models of fluid, hexatic and polymerized membranes. J. Phys.: Cond. Matter. 9, p. 8795
 P. B. S. Kumar, G. Gompper, and R. Lipowsky (2001) Budding dynamics of multicomponent membranes. Phys. Rev. Lett. 86, p. 3911
 P. B. S. Kumar and M. Rao (1998) Shape instabilities in the dynamics of a twocomponent fluid membrane. Phys. Rev. Lett. 80, p. 2489
 G. Gompper and D. M. Kroll (1997) Freezing flexible vesicles. Phys. Rev. Lett. 78, p. 2859
 G. Gompper and D. M. Kroll (1998) Membranes with fluctuating topology: Monte Carlo simulations. Phys. Rev. Lett. 81, p. 2284
 J. C. Shillcock and D. H. Boal (1996) Entropydriven instability and rupture of fluid membranes. Biophys. J. 71, p. 317
 H. Noguchi and G. Gompper (2005) Fluid vesicles with viscous membranes in shear flow. Phys. Rev. Lett. 93, p. 258102
 H. Noguchi and G. Gompper (2005) Shape transitions of fluid vesicles and red blood cells in capillary flows. PNAS 102, p. 14159
 H. Noguchi and G. Gompper (2006) Meshless membrane model based on the moving leastsquares method. Phys. Rev. E 73, p. 021903

G. T. Linke (2005) Dissertation Universität Potsdam.
URN: urn:nbn:de:kobv:517opus5835
URL: http://opus.kobv.de/ubp/volltexte/2005/583/.  G. T. Linke, R. Lipowsky, and T. Gruhn (2006) Osmotically induced passage of vesicles through narrow pores, Europhys. Lett. 74, p. 916
 E. Guyon, J.P. Hulin, L. Petit, and C. D. Mitescu (2001) Physical Hydrodynamics. Oxford University Press, Oxford
 E. B. Bagley, I. M. Cabott, and D. C. West (1958) Discontinuity in the flow curve of polyethylene. J. Appl. Phys. 29, p. 109
 T. C. B. McLeish and R. C. Ball (1986) A molecular approach to the spurt effect in polymer melt flow. J .Polym. Sci. 24, p. 1735
 M. E. Cates, T. C. B. McLeish, and G. Marrucci (1993) The rheology of entangled polymers at very high shear rates. Europhys. Lett. 21, p. 451
 D. C. Roux, J.F. Berret, G. Porte, E. PeuvrelDisdier, and P. Lindner (1995) Shearinduced orientations and textures of nematic wormlike micelles. Macromolecules 28, p. 1681
 J.F. Berret, G. Porte, and J.P. Decruppe (1996) Inhomogeneous shear rows of wormlike micelles: A master dynamic phase diagram. Phys. Rev. E 55, p. 1668
 M. M. Britton and P. T. Callaghan (1999) Shear banding instability in wormlike micellar solutions Eur. Phys. J. B 7, p. 237
 E. Fischer and P. T. Callaghan (2001) Shear banding and the isotropictonematic transition in wormlike micelles. Phys. Rev. E 64, p. 011501
 M. R. LopezGonzalex, W. M. Holmes, P. T. Callaghan, and P. J. Photinos (2004) Shear banding fluctuations and nematic order in wormlike micelles. Phys. Rev. Lett. 93, p. 268302
 P. D. Olmsted and P. M. Goldbart (1992) Isotropicnematic transition in shear flow: State selection, coexistence, phase transitions, and critical behavior. Phys. Rev. A 46, p. 4966
 G. Porte, J.F. Berret, and J. L. Harden (1997) Inhomogeneous flows of complex fluids: Mechanical instability versus nonequilibrium phase transition. J. de Physique II 7, p. 459
 V. Schmitt, C. M. Marques, and F. Lequeux (1995) Shearinduced phase separation of complex fluids  the role of flowconcentration coupling. Phys. Rev. E 52, p. 4009
 P. D. Olmsted and C. Y. D. Lu (1997) Coexistence and phase separation in sheared complex fluids. Phys. Rev. E 56, p. R55
 M. P. Lettinga and J. K. G. Dhont (2004) Nonequilibrium phase behaviour of rodlike viruses under shear flow. J. Phys.: Cond. Matter 16, p. S3929
 P. D. Olmsted and C. Y. D. Lu (1999) Phase separation of rigidrod suspensions in shear flow. Phys. Rev. E 60, p. 4397
 P. D. Olmsted (1999) Twostate shear diagrams for complex fluids in shear flow. Europhys. Lett. 48, p. 339
 S. M. Fielding and P. D. Olmsted (2003) Flow phase diagrams for concentrationcoupled shear banding. Europ. Phys. J. E 11, p. 65
 N. K. Ailawadi, B. J. Berne, and D. Forster (1971) Hydrodynamics and collective angularmomentum fluctuations in molecular fluids. Phys. Rev. A 3, p. 1462
 X. F. Yuan, M. P. Allen (1997) Nonlinear responses of the hardspheroid fluid under shear flow. Physica A 240, p. 145
 H. See, M. Doi, and R. Larson (1990) The effect of steady flow fields on the isotropicnematic phase transition of rigid rodlike polymers. J. Chem. Phys. 92, p. 792
 P. D. Olmsted and P. Goldbart (1990) Theory of the nonequilibrium phase transition for nematic liquid crystals under shear flow. Phys. Rev. A 41, p. 4578
 J. F. Berret, D. C. Roux, and G. Porte (1994) Isotropictonematic transition in wormlike micelles under shear. J. de Physique II 4, p. 1261
 J. F. Berret, D. C. Roux, G. Porte, and P. Lindner (1994) Shearinduced isotropictonematic phase transition in equilibrium polymers. Europhys. Lett. 25, p. 521
 E. Cappelaere, J.F. Berret, J. P. Decruppe, R. Cressely, and P. Lindner (1997) Rheology, birefringence, and smallangle neutron scattering in a charged micellar system: Evidence of a shearinduced phase transition. Phys. Rev. E 56, p. 1869
 J.F. Berret, D. C. Roux, and P. Lindner (1998) Structure and rheology of concentrated wormlike micelles at the shearinduced isotropictonematic transition. Eur. Phys. J. B 5, p. 67
 P. T. Mather, A .RomoUribe, C. D. Han, and S. S. Kim (1997) Rheooptical evidence of a flowinduced isotropicnematic transition in a thermotropic liquidcrystalline polymer. Macromolecules 30, p. 7977
 R. G. Larson and D. W. Mead (1993) The Ericksen number and Deborah number cascades in sheared polymeric nematics. Liquid Crystals 15, p. 151
 J. F. Berret, D. C. Roux, G. Porte, and P. Lindner (1995) Tumbling behavior of nematic wormlike micelles under shearflow Europhys. Lett. 32, p. 137
 A. V. Zakharov, A. A. Vakulenko, and J. Thoen (2003) Tumbling instability in a shearing nematic liquid crystal: Analysis of broadband dielectric results and theoretical treatment. J. Chem. Phys. 118, p. 4253
 S. Hess, M. Kröger (2004) Regular and chaotic orientational and rheological behaviour of liquid crystals. J. Phys.: Cond. Matter 16, p. S3835
 S. Sarman and D. J. Evans (1993) Statistical mechanics of viscous flow in nematic fluids. J. Chem. Phys. 99, p. 9021
 S. Tang, G. T. Evans, C. P. Mason, and M. P. Allen (1995) Shear viscosity for fluids of hard ellipsoids  A kinetictheory and moleculardynamics study. J. Chem. Phys. 102, p. 3794
 M. P. Allen and D. J. Tildesley (1989) Computer Simulation of Liquids. Oxford University Press, New York
 D. J. Evans and T. P. Morriss (1990) Statistical Mechanics of Nonequilibrium Fluids. Academic Press, San Diego
 S. S. Sarman, D. J. Evans, and P. T. Cummings (1992) Recent developments in nonNewtonian molecular dynamics. Physics Reports 305, p. 1
 F. Varnik and K. Binder (2002) Shear viscosity of a supercooled polymer melt via nonequilibrium molecular dynamics simulations. J. Chem. Phys. 117, p. 6336
 A. W. Lees and S. F. Edwards (1972) Computer study of transport processes under extreme conditions. J. Phys. C 5, p. 1921
 D. J. Evans and G. P. Morriss (1984) Nonlinearresponse theory for steady planar couetteflow. Phys. Rev. A 30, p. 1528
 B. J. Edwards, C. Baig, and D. J. Keffer (2005) An examination of the validity of nonequilibrium moleculardynamics simulation algorithms for arbitrary steadystate flows. J. Chem. Phys. 123, p. 114106
 F. Zhang, D. J. Searles, D. J. Evans, J. S. D Hansen, and D. J. Isbister (1999) Kinetic energy conserving integrators for Gaussian thermostatted SLLOD. J. Chem. Phys. 111, p. 18
 G. A. Pan, J. F. Ely, C. McCabe, and D. J. Isbister (2005) Operator splitting algorithm for isokinetic SLLOD molecular dynamics. J. Chem. Phys 122, p. 094114
 D. Baalss and S. Hess (1986) Nonequilibrium moleculardynamics studies on the anisotropic viscosity of perfectly aligned nematic liquid crystals. Phys. Rev. Lett. 57, p. 86
 S. Sarman (1995) Nonequilibrium molecular dynamics of liquidcrystal shearflow. J. Chem. Phys. 103, p. 10378
 S. Sarman (1997) Shear flow simulations of biaxial nematic liquid crystals. J. Chem. Phys. 107, p. 3144
 J. Liam McWhirter and G. N. Patey (2002) Nonequilibrium molecular dynamics simulations of a simple dipolar fluid under shear flow. J. Chem. Phys. 117, p. 2747
 J. Liam McWhirter and G. N. Patey (2002) Molecular dynamics simulations of a ferroelectric nematic liquid under shear flow. J. Chem. Phys. 117, p. 8551
 F. MüllerPlathe (1999) Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids. Phys. Rev. E 59, p. 4894
 D. Frenkel and B. Smit (2002) Understanding Molecular Simulations. Academic Press, San Diego
 D. J. Evans and S. Sarman (1993) Equivalence of thermostatted nonlinear responses. Phys. Rev. E 48, p. 65
 D. Ruelle (2000) A Remark on the Equivalence of Isokinetic and Isoenergetic Thermostats in the Thermodynamic Limit. J. Stat. Phys. 100, p. 757
 A. Kolb and B. Dünweg (1999) Optimized constant pressure stochastic dynamics. J. Chem. Phys. 111, p. 4453
 T. Soddemann, B. Dünweg, and K. Kremer (2003) Dissipative particle dynamics: A useful thermostat for equilibrium and nonequilibrium molecular dynamics simulations. Phys. Rev. E 68, p. 046702
 R. Zwanzig (1961) Memory effects in irreversible thermodynamics. Phys. Rev. 124, p. 983
 H. C. Öttinger (1998) General projection operator formalism for the dynamics and thermodynamics of complex fluids. Phys. Rev. E 57, p. 1416
 A. N. Gorban, I. V. Karlin, H. C. Öttinger, and L. L. Tatarinova (2001) Ehrenfests argument extended to a formalism of nonequilibrium thermodynamics. Phys. Rev. E 63, p. 066124
 B. Dünweg and W. Paul (1991) Brownian dynamics simulations without Gaussian random numbers. Int. J. Mod. Phys. C 2, p. 817
 G. Germano and F. Schmid (2003) Simulation of nematicisotropic phase coexistence in liquid crystals under shear. In Publication Series of the John von Neumann Institute for Computing, Vol. 20, p. 311
 G. Germano and F. Schmid (2005) Nematicisotropic interfaces under shear: A molecular dynamics simulation. J. Chem. Phys. 123, p. 214703
 M. P. Allen (2000) Molecular simulation and theory of the isotropicnematic interface. J. Chem. Phys. 112, p. 5447
 A. J. McDonald, M. P. Allen, and F. Schmid (2001) Surface tension of the isotropicnematic interface. Phys. Rev. E 63, p. 10701R
 N. Akino, F. Schmid, and M. P. Allen (2001) Moleculardynamics study of the nematicisotropic interface. Phys. Rev. E 63, p. 041706
 H. Guo, K. Kremer, and T. Soddemann (2002) Nonequilibrium molecular dynamics simulation of shearinduced alignment of amphiphilic model systems. Phys. Rev. E 66, p. 061503
 T. Soddemann, G. K. Auernhammer, H. Guo, B. Dünweg, and K. Kremer (2004) Shearinduced undulation of smecticA: Molecular dynamics simulations vs. analytical theory. Eur. Phys. J. E 13, p. 141
 V. K. Gupta, R. Krishnamoorti, J. A. Kornfield, and S. D. Smith (1995) Evolution of microstructure during shear alignment in a polystyrenepolyisoprene lamellar diblock copolymer. Macromolecules 28, p. 4464
 K. A. Koppi, M. Tirrell, F. S. Bates, K. Almdal, and R. H. Colby (1992) Lamellae orientation in dynamically sheared diblock copolymer melts. J. Physique II 2, p. 1941
 G. K. Auernhammer, H. R. Brand, and H. Pleiner (2000) The undulation instability in layered systems under shear flow  a simple model. Rheol. Acta 39, p. 215
 G. K. Auernhammer, H. R. Brand, and H. Pleiner (2002) Shearinduced instabilities in layered liquids. Phys. Rev. E 66, p. 061707
 A. D. Rey and M. M. Denn (2002) Dynamical phenomena in liquidcrystalline materials. Annu. Rev. Fluid Mech. 34, p. 233
 C. Denniston, E. Orlandini, and J. M. Yeomans (2000) Simulations of liquid crystal hydrodynamics in the isotropic and nematic phases. Europhys. Lett. 52, p. 481
 C. Denniston, D. Marenduzzo, E. Orlandini, and J. M. Yeomans (2004) Lattice Boltzmann algorithm for threedimensional liquidcrystal hydrodynamics. Phil. Trans. Royal Soc. London A 362, p. 1745
 A .N. Beris, B. J. Edwards, and M. Grmela (1990) Generalized constitutive equation for polymeric liquid crystals. 1. Model formulation using the Hamiltonian (Poisson bracket) formulation. J. NonNewton. Fluid Mechanics 35, p. 51
 D. Marenduzzo, E. Orlandini, and J. M. Yeomans (2003) Rheology of distorted nematic liquid crystals. Europhys. Lett. 64, p. 406
 D. Marenduzzo, E. Orlandini, and J. M. Yeomans (2004) Interplay between shear flow and elastic deformations in liquid crystals. J. Chem. Phys. 121, p. 582