1 Introduction

Coarse-grained models of complex fluids at equilibrium and unter shear

Friederike Schmid

Fakultät für Physik, Universität Bielefeld,

Universitätsstraße 25, D-33615 Bielefeld

E-mail: schmidphysik.uni-bielefeld.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 coarse-grained models at different levels of coarse-graining can provide insight into generic structures and basic dynamical processes at equilibrium and non-equilibrium.
In the first part of this lecture, some popular coarse-grained models for membranes and membrane systems are reviewed. Special focus is given to bead-spring models with different solvent representations, and to random-interface models. Selected examples of simulations at the molecular and the mesoscopic level are presented, and it is shown how simulations of molecular coarse-grained 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 non-equilibrium 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 non-Newtonian [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 small-scale degrees of freedom in the large-scale models, and to incorporate the details of the small-scale properties in the parameters of a new, simpler, model. This is called coarse-graining.

There are two aspects to coarse-graining. First, systematic coarse-graining 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, coarse-grained 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 coarse-grained 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 coarse-grained 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 coarse-grained approaches for membrane systems and liquid crystals.

2 Coarse-grained models for surfactant layers

In this section, we shall discuss coarse-grained simulation models and methods for amphiphilic membranes and membrane stacks. After a brief introduction into the phenomenology and some general remarks on coarse-graining, we will focus on two particular types of coarse-grained models: Particle-based bead-spring 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.

Figure 1: Artist’s view of a biomembrane. Source: NIST

This structure rests on a propensity of lipids to self-assemble 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].

Figure 2: Examples of membrane lipids. Left: Phosphatidylcholine (unsaturated); Right: Sphingomyelin.

By assembling into bilayers, the lipids can arrange themselves such that the head groups shield the tails from the water environment. Therefore, lipids in lipid-water 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.

Figure 3: Self-assembled bilayer structures.

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 lipid-water 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 one-component 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].

Figure 4: Some prominent membrane phases.

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:

  • Self-assembling 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 self-assembly: 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 coarse-graining. 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 coarse-grained models, up to mesoscopic models that are based on continuum theories for membranes.

In the following, we shall focus on two important examples: Coarse-grained 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 Coarse-grained molecular models

In coarse-grained 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 self-assembly), and the molecular flexibility (for the internal phase transitions).

Spring-bead models

The first coarse grained molecular models have been formulated on a lattice, relatively shortly after the introduction of the above-mentioned 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 self-assembly 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 off-lattice models.

Off-lattice molecular approaches to study self-assembling 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. Non-bonded beads interact via simple short-range pairwise potentials, e.g., a truncated Lennard-Jones potential. The parameters of the potentials are chosen such that pairs and pairs effectively repel each other. The Smit model reproduces self-assembly, micelle formation and membrane formation [33, 34, 35, 36, 37].

Most molecular off-lattice 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 non-bonded 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 two-dimensional 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 well-defined. 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 self-assembly is a non-trivial 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 solvent-free model. The crucial requirement seems to be that the interactions between hydrophobic beads have a smooth attractive part of relatively long range [50].

Figure 5: Schematic picture of bead-spring models for self-assembling membranes. (a) Explicit solvent model (b) Sandwich model (c) Solvent free model (d) Phantom solvent model.

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 self-assembled. Thus the self-assembly is to a large extent driven entropically.

In Monte Carlo simulations, this model is just as efficient as solvent-free 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 self-assembly. 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 Smit-type 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 self-assembled membrane, and discussed the advantages and disadvantages. The different models are sketched schematically in Fig. 5.

We shall now illustrate the potential of molecular coarse-grained models with two examples.

First application example: The ripple phase

Our first example shows how coarse-grained 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 freeze-etch micrographs, AFM pictures, and X-ray 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 X-ray 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.

Figure 6: Electron density maps of (top) an asymmetric rippled state (DMPC at 18.2 degrees) and (bottom) a symmetric rippled state (DPPC at 39.2 ). Reprinted with permission from Sengupta, Raghunathan, Katsaras [53] (rescaled). Copyright 2003 by the American Physical Society.

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

Figure 7: Asymmetric ripple in a Lecithin Bilayer, as observed in an atomistic simulation. From de Vries et al [54]. Copyright 2005 National Academy of Sciences, U.S.A.

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 coarse-grained models.

Rippled phases had also been observed in simulations of coarse-grained 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 bead-spring 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 Lennard-Jones 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.

Figure 8: Ripples from simulations of a coarse-grained model. (a) Asymmetric ripple (b) Symmetric ripple.

Membranes were found to self-assemble 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 coarse-grained 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 coarse-grained membrane models can be used to test and verify mesoscopic models, and how they can bridge between coarse-graining 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


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 coarse-grained molecular simulations (C. Loison et al [60]) of a binary mixture of amphiphiles and solvent, within a Smit-type 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:

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

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

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

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

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


Figure 9: Left: Snapshot of a bilayer stack (30720 amphiphiles and 30720 solvent beads). The “hydrophobic” tail beads are dark, the “hydrophilic” head beads and the solvent beads are light; Right: Snapshot of a single membrane position .

Figure 10: Ratios of transmembrane structure factors and vs. in-plane wavevector in units of . The solid lines correspond to a theoretical fit to Eq. (2) with one (common) fit parameter . After Ref. [60].

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]


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.

Figure 11: Snapshot of a single bilayer (top view. Only hydrophobic beads are shown).

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]


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.


Figure 12: Distribution of pore contours in a semi-logarithmic plot. Left: Raw data, Right: Divided by the degeneracy function . After Ref. [67].

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


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 two-dimensional self-avoiding 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 self-avoiding 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


Figure 13 shows that this is indeed the case.

Figure 13: Pore area vs. contour length (arbitrary units). After Ref. [67].

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 coarse-grained 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, coarse-grained 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 random-interface theories.

Random interface models

In random-interface models, membranes are represented by two-dimensional 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]


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 Gauss-Bonnet theorem states that


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 self-avoiding, 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 self-avoiding surfaces, and second a way to discretize the free energy functional (6) on those surfaces.

A widely used method to generate self-avoiding 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 (tethered-bead model). The tethers define the neighbor relations on the surface. In order to enforce self-avoidance, 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.

Figure 14: Connectivity update in the dynamic triangulation algorithm

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 Euler-characteristic, 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 non-trivial, 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


This completes the tethered-bead 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 tethered-bead 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 self-assemble into a two-dimensional membrane without being connected by tethers. The membrane is coarse-grained as a curved surface as in tethered-bead 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 tethered-bead 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.

Figure 15: Osmotically driven vesicle translocation through a pore. Configuration snapshots. From G. T. Linke [88].

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


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, pre-determined 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.

Figure 16: Free energy landscape on crossing the pore for different osmolarities starting from (open circles) to (plus) in steps of 100. All units are arbitrary. The osmolarity is kept fixed. See text for more explanation. From G. Linke et al[88, 89].

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 particle-based model would have been much more expensive, and both the setup and the analysis would have been more difficult: One has no well-defined 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 coarse-grained models for membrane systems: Bead-spring models for simulations on the coarse-grained molecular level, and random interface models for the coarse-grained 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 particle-based 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 coarse-grained simulations of complex fluids under shear. As in the previous section, we will move from the coarse-grained 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 off-diagonal components give the shear stress and correspond to tangential forces.

A fluid element with the velocity is thus accelerated according to


where is the local mass density. The stress tensor is often divided in three parts.


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 so-called Newtonian fluids, which have no elasticity () and a linear, instantaneous relation between the viscous stress tensor and the velocity gradient,


In the case of an incompressible, isotropic, Newtonian fluid, one has only one independent material parameter, the kinematic viscosity , and the constitutive equation reads


where the value of the local pressure is determined by the requirement that the incompressibility condition


is always fulfilled. Equations (10) in combination with (13) are the famous Navier-stokes equations (without external forces).

Complex fluids have internal degrees of freedom, and large characteristic length and time scales. Therefore, they become non-Newtonian 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 stress-strain 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 stress-strain 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).


Figure 17: Sketches of nonlinear stress-strain flow curves Left: Shear thickening and shear thinning; Right: Nonmonotonic stress-strain relation leading to phase separation. The horizontal dotted lines connect the two coexisting phases for three selected global strain rates.

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.

Figure 18: Examples of liquid crystalline phases. In the nematic phase (N), the particles are oriented, but translationally disordered. In the smectic phases, they form layers in one direction, but remain fluid within the layers. In the smectic B and smectic F phase, the structure within the layers is not entirely disordered, but has a type of order called hexatic. See Refs. [15, 16] for explanation. The isotropic phase (I) is entirely disordered (regular isotropic liquid).

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


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]


Figure 19 illustrates the corresponding fundamental distortions, the splay (), twist (), and bend mode ().

Figure 19: Elastic modes in nematic liquid crystals.

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


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 shear-induced 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 stress-strain 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


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]


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 so-called Leslie coefficients. They are linked by one relation [16],


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,


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 ,


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 Leslie-Ericksen equations of nemato-hydrodynamics [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 Ericksen-Leslie 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 stress-strain relationship, which eventually lead to mechanical phase separation [99].

Figure 20 shows an experimental example of a nonequilibrium phase diagram for a system of rodlike virusses [103].

Figure 20: Nonequilibrium phase diagram of a lyotropic liquid crystal under shear as a function of strain rate and particle density . ( is the density of the equilibrium isotropic/nematic transition). The system is a mixture of rodlike virusses (fd virus) and polymers (dextran). From Ref. [103].

3.2 Simulating shear on the particle level: NEMD

In the linear response regime, i.e., in the low-shear 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 Green-Kubo relations [120, 121].

However, the non-Newtonian character of most complex fluids and their resulting unique properties manifest themselves only beyond the linear response regime. To study these, non-equilibrium simulations are necessary. We will now briefly review non-equilibrium 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.

Figure 21: Lees-Edwards boundary conditions.

An alternative straightforward way to generate planar Couette flow is to use moving periodic boundary conditions as illustrated in Fig. 21 (Lees-Edwards 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,


In the -direction, the replicated particles acquire an additional velocity and an offset in the -direction,


Here are the dimensions of the simulation box, and the unit vector in the -direction.

Lees-Edwards 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


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 operator-split 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üller-Plathe [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 zigzag-shaped shear profile (Fig. 22).

Figure 22: Schematic sketch of the Müller-Plathe method for imposing shear. In periodical intervals, the particle in the middle slab that moves fastest in the -direction and the particle in the top slab that moves fastest in the ()-direction exchange their momentum components . After Ref. [136].

The Müller-Plathe 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 W-shaped 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üller-Plathe 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 pre-imposed profile ( in our earlier examples). Thermostats that use the pre-imposed flow profile are called “profile-biased thermostats”. At high strain rates, and in inhomogeneous systems, the use of profile-biased 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 “profile-unbiased 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 so-called 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


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 well-known 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 well-established 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 coarse-grained model, but can still absorb and carry heat. From the theory of coarse-graining, it is well-known 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 ,


where fulfills


i.e., the noise is delta-correlated and Gaussian distributed with mean zero and variance (fluctuation-dissipation 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 time-consuming. In fact, even in coarse-grained molecular models, simulated strain rates are usually much higher than experimentally accessible strain rates.

We now give two examples of recent large-scale NEMD studies of inhomogeneous complex systems under shear.

First application example: Nematic-Isotropic Interfaces under shear

Our first example is a study of the behavior of a Nematic-Isotropic 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].


Figure 23: Snapshot of a sheared nematic-isotropic interface (115200 particles) at strain rate . Also shown is the coordinate system used throughout this section. From Ref. [146].

Shear was imposed with Lees-Edwards 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 flow-aligned setup, the director points in the direction of flow (the direction); in the log-rolling setup, it points in the direction of “vorticity” (the direction). Only the flow-aligned 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 coarse-grained profiles by convoluting the profile with a symmetrical box-like coarse-graining function of width . The coarse-graining width must be chosen such that the coarse-grained 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.

Figure 24: Illustration of block analysis to obtain local profiles.

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.

Figure 25: Interfacial order parameter profile at strain rate (flow aligned setup). Inset shows the difference between profiles with and without shear. From Ref. [147].

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.


Figure 26: Velocity profile in the flow direction (left), and in the vorticity direction (right), at strain rate . From Ref. [147].

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 .

Figure 27: Phase diagram of the flow-aligned system. The dashed tie lines connect coexisting states. I denotes the isotropic region, and N the nematic region. From Ref. [147].

Second application example: Shear-induced phenomena in surfactant layers

As a second example, we discuss the effect of shear on lamellar stacks of surfactants. This has been studied in coarse-grained 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üller-Plathe 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.

Figure 28: Reorientation of surfactant layers under shear from an initial transverse state at low strain rate, . Reprinted with permission from Ref. [151]. Copyright 2002 by the American Physical Society.

Figure 29: Reorientation of surfactant layers under shear from an initial transverse state at high strain rate, . Reprinted with permission from Ref. [151]. Copyright 2002 by the American Physical Society.


Figure 30: Shear induced undulations. Left: configuration snapshot at strain rate . Right: Undulation amplitude as a function of strain rate. Undulations set in at a strain rate of . From Ref. [152].

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 well-defined strain-rate.

Our two examples show that coarse-grained 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 Leslie-Ericksen 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 Leslie-Ericksen 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 liquid-crystal hydrodynamics, which allows to simulate liquid crystal flow [158, 159]. The starting point is the Beris-Edwards 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 Lattice-Boltzmann 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].

Figure 31: Snapshots from a Lattice-Boltzmann simulation of a sheared hybrid aligned nematic cell. White regions are ordered, and black regions are disordered. The numbers give time in milliseconds. From Ref. [161].

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 (no-slip 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 flow-aligned, 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 coarse-graining have been developed. In this section, we have presented and discussed several variants of non-equilibrium 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 coarse-grained 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 coarse-grained 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.


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, Siewert-Jan 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.


  1. M. Daoud and C. E. Williams, Eds. (1995) Soft Matter Physics. Springer-Verlag, Berlin
  2. R. G. Larson (1999) The Structure and Rheology of Complex Fluids. Oxford University Press, New York
  3. J. N. Israelachvili (1992) Intermolecular and Surface Forces. Academic Press, London
  4. M. Kröger (2003) Simple models for complex nonequilibrium fluids. Physics reports 390, p. 453
  5. P. J. Flory (1969) Statistical Mechanics of Chain Molecules. Interscience Publishers, New York
  6. P.-G. de Gennes (1979) Scaling Concepts in Polymer Physics. Cornell University Press, Ithaca
  7. M. Doi (1992) Introduction to Polymer Physics. Clarendon Press, Oxford
  8. M. Doi and S. Edwards (1986) The Theory of Polymer Dynamics. Clarendon Press, Oxford
  9. W. B. Russel, D. A. Saville, and W. R. Schowalter (1989) Colloidal Dispersions. Cambridge University Press, Cambridge
  10. R. J. Hunter (1989) Foundations of Colloid Science. Clarendon Press, Oxford
  11. M. Borowko, Ed. (2000) Computational Methods in Surface and Colloid Science. Marcel Dekker Inc., New York
  12. S. A. Safran (1994) Statistical Thermodynamics of Surfaces, Interfaces, and Membranes. Addison-Wesley, Reading, MA
  13. G. Gompper and M. Schick (1994) Self-Assembling Amphiphilic Systems, in Phase Transitions and Critical Phenomena, C. Domb and J. Lebowitz Eds., Academic Press, London, 1994 16, p. 1.
  14. P. G. de Gennes and J. Prost (1995) The Physics of Liquid Crystals. Oxford University Press, Oxford
  15. S. Chandrasekhar (1992) Liquid Crystals. Cambridge University Press, Cambridge
  16. P. M. Chaikin and T. C. Lubensky (1995) Principles of Condensed Matter Physics. Cambridge University Press, Cambridge
  17. F. Müller-Plathe (2002) Coarse-graining in polymer simulation: From the atomistic to the mesoscopic scale and back. ChemPhysChem 3, p. 754
  18. 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
  19. R. B. Gennis (1989) Biomembranes. Springer Verlag, New York
  20. 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
  21. O. G. Mouritsen (2005) Life – As a Matter of Fat. Springer, Berlin Heidelberg
  22. K. V. Schubert and E. W. Kaler (1996) Nonionic microemulsions. Ber. Bunseng. Phys. Chemie 100, p. 190
  23. S. Doniach (1978) Thermodynamic fluctuations in phospholipid bilayers. J. Chem. Phys. 68, p. 4912
  24. D. A. Pink, T. J. Green, and D. Chapman (1980) Raman-scattering in bilayers of saturated phosphatidylcholines - experiment and theory. Biochemistry 19, p. 349
  25. A. Caillé, D. A. Pink, F. de Verteuil, and M. Zuckermann (1980) Theoretical models for quasi-two-dimensional mesomorphic monolayers and membrane bilayers. Can. J. Physique 58, p. 581
  26. 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.
  27. 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
  28. R. G. Larson, L. E. Scriven, and H. T. Davis (1985) Monte-Carlo simulation of model amphiphilic oil-water systems. J. Chem. Phys. 83, p. 2411
  29. 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.
  30. R. G. Larson (1996) Monte Carlo simulations of the phase behavior of surfactant solutions. J. Physique II 6, p. 1441
  31. B. Smit, A. G. Schlijper, L. A. M. Rupert, and N. M. van Os (1990) Effects of chain-length of surfactants on the interfacial tension - molecular-dynamics simulations and experiments. J. Phys. Chem. 94, p. 6933
  32. B. Smit, K. Esselink, P. A. J. Hilbers, N. M. van Os, L. A. M. Rupert, and I. Szleifer (1993) Computer simulations of surfactant self-assembly. Langmuir 9, p. 9
  33. S. Karaborni, K. Esselink, P. A. J. Hilbers, B. Smit, J. Karthäuser, N. M. van Os, and R. Zana, Simulating the self-assembly of gemini (dimeric) surfactants. Science 266, p. 254
  34. B. J. Palmer and J. Liu (1996) Simulations of micelle self-assembly in surfactant solutions. Langmuir 12, p. 746
  35. B. J. Palmer and J. Liu (1996) Effects of solute-surfactant interactions on micelle formation in surfactant solutions Langmuir 12, p. 6015
  36. R. Goetz and R. Lipowsky (1998) Computer simulations of bilayer membranes: Self-assembly and interfacial tension. J. Chem. Phys. 108, p. 7397
  37. 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
  38. M. J. Stevens (2004) Coarse-grained simulations of lipid bilayers . J. Chem. Phys. 121, p. 11942
  39. M. Kranenburg, J. P. Nicolas, and B. Smit (2004) Comparison of mesoscopic phospholipid-water models. Phys. Chem. Chem. Phys. 6, p. 4142
  40. A. F. Jakobsen, O. G. Mouritsen, and G. Besold (2005) Artifacts in dynamical simulations of coarse-grained model lipid bilayers. J. Chem. Phys. 122, p. 204901
  41. D. Harries and A. Ben-Shaul (1997) Conformational chain statistics in a model lipid bilayer: Comparison between mean field and Monte Carlo calculations. J. Chem. Phys. 106, p. 1609
  42. A. Baumgärtner (1995) Asymmetric partitioning of a polymer into a curved membrane. J. Chem. Phys. 103, p. 10669
  43. A. Baumgärtner (1996) Insertion and hairpin formation of membrane proteins: A Monte Carlo study. Biophys. J. 71, p. 1248
  44. T. Sintes and A. Baumgärtner (1997) Short-range attractions between two colloids in a lipid monolayer. Biophys. J. 73, p. 2251
  45. O. Lenz and F. Schmid (2004) A simple computer model for liquid lipid bilayers. J. Mol. Liquids 117, p. 147
  46. R. Goetz, G. Gompper, and R. Lipowsky, Mobility and elasticity of self-assembled membranes. Phys. Rev. Lett. 82, p. 221
  47. H. Noguchi and M. Takasu (2001) Self-assembly of amphiphiles into vesicles: A Brownian dynamics simulation. Phys. Rev. E 64, p. 041913
  48. O. Farago (2003) ”Water-free” computer model for fluid bilayer membranes. J. Chem. Phys. 119, p. 596
  49. I. R. Cooke, K. Kremer, and M. Deserno (2005) Tunable generic model for fluid bilayer membranes. Phys. Rev. E 72, p. 011506
  50. H. Noguchi and M. Takasu (2001) Fusion pathways of vesicles: A Brownian dynamics simulation. J. Chem. Phys. 115, p. 9547
  51. H. Noguchi and M. Takasu (2002) Structural changes of pulled vesicles: A Brownian dynamics simulation. Phys. Rev. E 65, p. 051907
  52. K. Sengupta, V. A. Raghunathan, and J. Katsaras (2003), Structure of the ripple phase of phospholipid multibilayers. Phys. Rev. E 68, p. 031710
  53. 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
  54. M. Kranenburg, C. Laforge, and B. Smit (2004) Mesoscopic simulations of phase transitions in lipid bilayers. Phys. Chem. Chem. Phys. 6, p. 4531
  55. M. Kranenburg and B. Smit (2005) Phase behavior of model lipid bilayers. J. Phys. Chem. B 109, p. 6553
  56. O. Lenz and F. Schmid (2007) Structure of symmetric and asymmetric ‘ripple’ phases in lipid bilayers. Phys. Rev. Lett. 98, p. 058104
  57. 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
  58. 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
  59. C. Loison, M. Mareschal, K. Kremer, and F. Schmid (2003) Thermal fluctuations in a lamellar phase of a binary amphiphile-solvent mixture: A molecular dynamics study. J. Chem. Phys. 119, p. 13138
  60. T. Soddemann, B. Dünweg, and K. Kremer (2001) A generic computer model for amphiphilic systems. Eur. Phys. J. E 6, p. 409
  61. M. Müller and M. Schick (1996) New mechanism of membrane fusion. J. Chem. Phys. 116, p. 2342
  62. S.-J. Marrink, F. Jähning, and H. Berendsen (1996) Proton transport across transient single-file water pores in a lipid membrane studied by molecular dynamics simulations. Biophys. J. 71, p. 632
  63. D. Zahn and J. Brickmann (2002) Molecular Dynamics Study of Water Pores in a Phospholipid Bilayer. Chem. Phys. Lett. 352, p. 441
  64. 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
  65. Z. J. Wang and D. Frenkel (2005) Pore nucleation in mechanically stretched bilayer membranes. J. Chem. Phys. 123, p. 154701
  66. C. Loison, M. Mareschal, and F. Schmid (2004) Pores in bilayer membranes of amphiphilic molecules: Coarse-Grained Molecular Dynamics Simulations Compared with Simple Mesoscopic Models. J. Chem. Phys. 121, p. 1890
  67. C. Loison, M. Mareschal, and F. Schmid (2005) Fluctuations and defects in lamellar stacks of amphiphilic bilayers. Comp. Phys. Comm. 169, p. 99
  68. J. D. Lister (1975) Stability of lipid bilayers and red blood-cell membranes. Physics Lett. A A 53, p. 193
  69. E. W. Weisstein (2003) CRC Concise encyclopaedia of mathematics. Chapman & Hall CRC, http://mathworld.wolfram.com
  70. W. Helfrich (1973) Elastic properties of lipid bilayers - theory and possible experiments. Z. Naturforschung C 28, p. 693
  71. E. Evans (1974) Bending resistance and chemically induced moments in membrane bilayers. Biophys. J. 14, p. 923
  72. U. Seifert (1997) Configurations of fluid membranes and vesicles. Adv. Phys. 46, p. 13
  73. 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
  74. A. Billoire and F. David (1986) Scaling properties of randomly triangulated planar random surfaces - a numerical study. Nucl. Phys. B 275, p. 617
  75. Y. Kantor, M. Kardar, and D. R. Nelson (1986) Statistical mechanics of tethered surfaces. Phys. Rev. Lett. 57, p. 791
  76. J.-S Ho, A. Baumgärtner (1990) Simulations of fluid self-avoiding membranes. Europhys. Lett. 12, p. 295
  77. D. M. Kroll and G. Gompper (1992) The conformations of fluid membranes - Monte-Carlo simulations. Science 255, p. 968
  78. G. Gompper and D. M. Kroll (1997) Network models of fluid, hexatic and polymerized membranes. J. Phys.: Cond. Matter. 9, p. 8795
  79. P. B. S. Kumar, G. Gompper, and R. Lipowsky (2001) Budding dynamics of multicomponent membranes. Phys. Rev. Lett. 86, p. 3911
  80. P. B. S. Kumar and M. Rao (1998) Shape instabilities in the dynamics of a two-component fluid membrane. Phys. Rev. Lett. 80, p. 2489
  81. G. Gompper and D. M. Kroll (1997) Freezing flexible vesicles. Phys. Rev. Lett. 78, p. 2859
  82. G. Gompper and D. M. Kroll (1998) Membranes with fluctuating topology: Monte Carlo simulations. Phys. Rev. Lett. 81, p. 2284
  83. J. C. Shillcock and D. H. Boal (1996) Entropy-driven instability and rupture of fluid membranes. Biophys. J. 71, p. 317
  84. H. Noguchi and G. Gompper (2005) Fluid vesicles with viscous membranes in shear flow. Phys. Rev. Lett. 93, p. 258102
  85. H. Noguchi and G. Gompper (2005) Shape transitions of fluid vesicles and red blood cells in capillary flows. PNAS 102, p. 14159
  86. H. Noguchi and G. Gompper (2006) Meshless membrane model based on the moving least-squares method. Phys. Rev. E 73, p. 021903
  87. G. T. Linke (2005) Dissertation Universität Potsdam.
    URN: urn:nbn:de:kobv:517-opus-5835
    URL: http://opus.kobv.de/ubp/volltexte/2005/583/.
  88. G. T. Linke, R. Lipowsky, and T. Gruhn (2006) Osmotically induced passage of vesicles through narrow pores, Europhys. Lett. 74, p. 916
  89. E. Guyon, J.-P. Hulin, L. Petit, and C. D. Mitescu (2001) Physical Hydrodynamics. Oxford University Press, Oxford
  90. E. B. Bagley, I. M. Cabott, and D. C. West (1958) Discontinuity in the flow curve of polyethylene. J. Appl. Phys. 29, p. 109
  91. 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
  92. 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
  93. D. C. Roux, J.-F. Berret, G. Porte, E. Peuvrel-Disdier, and P. Lindner (1995) Shear-induced orientations and textures of nematic wormlike micelles. Macromolecules 28, p. 1681
  94. 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
  95. M. M. Britton and P. T. Callaghan (1999) Shear banding instability in wormlike micellar solutions Eur. Phys. J. B 7, p. 237
  96. E. Fischer and P. T. Callaghan (2001) Shear banding and the isotropic-to-nematic transition in wormlike micelles. Phys. Rev. E 64, p. 011501
  97. M. R. Lopez-Gonzalex, 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
  98. P. D. Olmsted and P. M. Goldbart (1992) Isotropic-nematic transition in shear flow: State selection, coexistence, phase transitions, and critical behavior. Phys. Rev. A 46, p. 4966
  99. G. Porte, J.-F. Berret, and J. L. Harden (1997) Inhomogeneous flows of complex fluids: Mechanical instability versus non-equilibrium phase transition. J. de Physique II 7, p. 459
  100. V. Schmitt, C. M. Marques, and F. Lequeux (1995) Shear-induced phase separation of complex fluids - the role of flow-concentration coupling. Phys. Rev. E 52, p. 4009
  101. P. D. Olmsted and C. Y. D. Lu (1997) Coexistence and phase separation in sheared complex fluids. Phys. Rev. E 56, p. R55
  102. M. P. Lettinga and J. K. G. Dhont (2004) Non-equilibrium phase behaviour of rod-like viruses under shear flow. J. Phys.: Cond. Matter 16, p. S3929
  103. P. D. Olmsted and C. Y. D. Lu (1999) Phase separation of rigid-rod suspensions in shear flow. Phys. Rev. E 60, p. 4397
  104. P. D. Olmsted (1999) Two-state shear diagrams for complex fluids in shear flow. Europhys. Lett. 48, p. 339
  105. S. M. Fielding and P. D. Olmsted (2003) Flow phase diagrams for concentration-coupled shear banding. Europ. Phys. J. E 11, p. 65
  106. N. K. Ailawadi, B. J. Berne, and D. Forster (1971) Hydrodynamics and collective angular-momentum fluctuations in molecular fluids. Phys. Rev. A 3, p. 1462
  107. X. F. Yuan, M. P. Allen (1997) Non-linear responses of the hard-spheroid fluid under shear flow. Physica A 240, p. 145
  108. H. See, M. Doi, and R. Larson (1990) The effect of steady flow fields on the isotropic-nematic phase transition of rigid rod-like polymers. J. Chem. Phys. 92, p. 792
  109. 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
  110. J. F. Berret, D. C. Roux, and G. Porte (1994) Isotropic-to-nematic transition in wormlike micelles under shear. J. de Physique II 4, p. 1261
  111. J. F. Berret, D. C. Roux, G. Porte, and P. Lindner (1994) Shear-induced isotropic-to-nematic phase transition in equilibrium polymers. Europhys. Lett. 25, p. 521
  112. E. Cappelaere, J.-F. Berret, J. P. Decruppe, R. Cressely, and P. Lindner (1997) Rheology, birefringence, and small-angle neutron scattering in a charged micellar system: Evidence of a shear-induced phase transition. Phys. Rev. E 56, p. 1869
  113. J.-F. Berret, D. C. Roux, and P. Lindner (1998) Structure and rheology of concentrated wormlike micelles at the shear-induced isotropic-to-nematic transition. Eur. Phys. J. B 5, p. 67
  114. P. T. Mather, A .Romo-Uribe, C. D. Han, and S. S. Kim (1997) Rheo-optical evidence of a flow-induced isotropic-nematic transition in a thermotropic liquid-crystalline polymer. Macromolecules 30, p. 7977
  115. R. G. Larson and D. W. Mead (1993) The Ericksen number and Deborah number cascades in sheared polymeric nematics. Liquid Crystals 15, p. 151
  116. J. F. Berret, D. C. Roux, G. Porte, and P. Lindner (1995) Tumbling behavior of nematic worm-like micelles under shear-flow Europhys. Lett. 32, p. 137
  117. 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
  118. S. Hess, M. Kröger (2004) Regular and chaotic orientational and rheological behaviour of liquid crystals. J. Phys.: Cond. Matter 16, p. S3835
  119. S. Sarman and D. J. Evans (1993) Statistical mechanics of viscous flow in nematic fluids. J. Chem. Phys. 99, p. 9021
  120. S. Tang, G. T. Evans, C. P. Mason, and M. P. Allen (1995) Shear viscosity for fluids of hard ellipsoids - A kinetic-theory and molecular-dynamics study. J. Chem. Phys. 102, p. 3794
  121. M. P. Allen and D. J. Tildesley (1989) Computer Simulation of Liquids. Oxford University Press, New York
  122. D. J. Evans and T. P. Morriss (1990) Statistical Mechanics of Nonequilibrium Fluids. Academic Press, San Diego
  123. S. S. Sarman, D. J. Evans, and P. T. Cummings (1992) Recent developments in non-Newtonian molecular dynamics. Physics Reports 305, p. 1
  124. F. Varnik and K. Binder (2002) Shear viscosity of a supercooled polymer melt via nonequilibrium molecular dynamics simulations. J. Chem. Phys. 117, p. 6336
  125. A. W. Lees and S. F. Edwards (1972) Computer study of transport processes under extreme conditions. J. Phys. C 5, p. 1921
  126. D. J. Evans and G. P. Morriss (1984) Nonlinear-response theory for steady planar couette-flow. Phys. Rev. A 30, p. 1528
  127. B. J. Edwards, C. Baig, and D. J. Keffer (2005) An examination of the validity of nonequilibrium molecular-dynamics simulation algorithms for arbitrary steady-state flows. J. Chem. Phys. 123, p. 114106
  128. 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
  129. 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
  130. D. Baalss and S. Hess (1986) Nonequilibrium molecular-dynamics studies on the anisotropic viscosity of perfectly aligned nematic liquid crystals. Phys. Rev. Lett. 57, p. 86
  131. S. Sarman (1995) Nonequilibrium molecular dynamics of liquid-crystal shear-flow. J. Chem. Phys. 103, p. 10378
  132. S. Sarman (1997) Shear flow simulations of biaxial nematic liquid crystals. J. Chem. Phys. 107, p. 3144
  133. 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
  134. 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
  135. F. Müller-Plathe (1999) Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids. Phys. Rev. E 59, p. 4894
  136. D. Frenkel and B. Smit (2002) Understanding Molecular Simulations. Academic Press, San Diego
  137. D. J. Evans and S. Sarman (1993) Equivalence of thermostatted nonlinear responses. Phys. Rev. E 48, p. 65
  138. D. Ruelle (2000) A Remark on the Equivalence of Isokinetic and Isoenergetic Thermostats in the Thermodynamic Limit. J. Stat. Phys. 100, p. 757
  139. A. Kolb and B. Dünweg (1999) Optimized constant pressure stochastic dynamics. J. Chem. Phys. 111, p. 4453
  140. 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
  141. R. Zwanzig (1961) Memory effects in irreversible thermodynamics. Phys. Rev. 124, p. 983
  142. H. C. Öttinger (1998) General projection operator formalism for the dynamics and thermodynamics of complex fluids. Phys. Rev. E 57, p. 1416
  143. 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
  144. B. Dünweg and W. Paul (1991) Brownian dynamics simulations without Gaussian random numbers. Int. J. Mod. Phys. C 2, p. 817
  145. G. Germano and F. Schmid (2003) Simulation of nematic-isotropic phase coexistence in liquid crystals under shear. In Publication Series of the John von Neumann Institute for Computing, Vol. 20, p. 311
  146. G. Germano and F. Schmid (2005) Nematic-isotropic interfaces under shear: A molecular dynamics simulation. J. Chem. Phys. 123, p. 214703
  147. M. P. Allen (2000) Molecular simulation and theory of the isotropic-nematic interface. J. Chem. Phys. 112, p. 5447
  148. A. J. McDonald, M. P. Allen, and F. Schmid (2001) Surface tension of the isotropic-nematic interface. Phys. Rev. E 63, p. 10701R
  149. N. Akino, F. Schmid, and M. P. Allen (2001) Molecular-dynamics study of the nematic-isotropic interface. Phys. Rev. E 63, p. 041706
  150. H. Guo, K. Kremer, and T. Soddemann (2002) Nonequilibrium molecular dynamics simulation of shear-induced alignment of amphiphilic model systems. Phys. Rev. E 66, p. 061503
  151. T. Soddemann, G. K. Auernhammer, H. Guo, B. Dünweg, and K. Kremer (2004) Shear-induced undulation of smectic-A: Molecular dynamics simulations vs. analytical theory. Eur. Phys. J. E 13, p. 141
  152. V. K. Gupta, R. Krishnamoorti, J. A. Kornfield, and S. D. Smith (1995) Evolution of microstructure during shear alignment in a polystyrene-polyisoprene lamellar diblock copolymer. Macromolecules 28, p. 4464
  153. 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
  154. 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
  155. G. K. Auernhammer, H. R. Brand, and H. Pleiner (2002) Shear-induced instabilities in layered liquids. Phys. Rev. E 66, p. 061707
  156. A. D. Rey and M. M. Denn (2002) Dynamical phenomena in liquid-crystalline materials. Annu. Rev. Fluid Mech. 34, p. 233
  157. 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
  158. C. Denniston, D. Marenduzzo, E. Orlandini, and J. M. Yeomans (2004) Lattice Boltzmann algorithm for three-dimensional liquid-crystal hydrodynamics. Phil. Trans. Royal Soc. London A 362, p. 1745
  159. 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. Non-Newton. Fluid Mechanics 35, p. 51
  160. D. Marenduzzo, E. Orlandini, and J. M. Yeomans (2003) Rheology of distorted nematic liquid crystals. Europhys. Lett. 64, p. 406
  161. 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
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.