Self-Assembly of Ternary Cubic, Hexagonal, and Lamellar

Feb 21, 2008 - We use a kinetic lattice-Boltzmann method to simulate the self-assembly of the cubic primitive (P), diamond (D), and gyroid (G) mesopha...
0 downloads 0 Views 764KB Size
Download PDF
2950

J. Phys. Chem. B 2008, 112, 2950-2957

Self-Assembly of Ternary Cubic, Hexagonal, and Lamellar Mesophases Using the Lattice-Boltzmann Kinetic Method R. S. Saksena† and P. V. Coveney* Centre for Computational Science, Department of Chemistry, UniVersity College London, 20 Gordon Street, London, WC1H 0AJ United Kingdom ReceiVed: April 24, 2007; In Final Form: NoVember 28, 2007

We use a kinetic lattice-Boltzmann method to simulate the self-assembly of the cubic primitive (P), diamond (D), and gyroid (G) mesophases from an initial quench composed of oil, water, and amphiphilic particles. Here, we also report the self-assembly of the noncubic hexagonal phase and two lamellar phases, one with periodic convolutions. The periodic mesophase structures are emergent from the underlying conservation laws and quasi-molecular interactions of the lattice-Boltzmann model. We locate regions of the model’s parameter space where the sequence of appearance of mesophases lamellar f primitive f hexagonal is in agreement with pressure jump experiments and the sequence cubic f lamellar is in agreement with compositional variations reported in the literature. The ability of our lattice-Boltzmann model to simulate self-assembly of cubic and noncubic phases in a unified and consistent manner opens the way for further investigations into the transition pathways and kinetics of the phase transitions between these states as well as of the rheology of these phases.

1. Introduction Mixtures of oil, water, and amphiphiles are known to selfassemble into a range of structures which are periodic in one, two, and three dimensions. In ternary amphiphilic mixtures, the amphiphiles organize themselves at the interfaces with their hydrophobic and hydrophilic ends oriented toward oil and water domains, respectively. Among the periodic surfaces observed in ternary amphiphilic mixtures, of particular interest in biotechnological applications are the triply periodic cubic bicontinuous phases defined by a minimal surface of zero curvature and with two interpenetrating fluid domains. Thermodynamically stable and long-lived metastable cubic phases have been observed for a range of surfactant solutions containing anionic and cationic soaps, nonionic and zwitterionic surfactants, as well as in solutions containing lipid molecules, the building blocks of biological membranes, and mixtures of lipids.1,2 Threedimensional cubic phases are known to be more stable than the one-dimensional lamellar and two-dimensional hexagonal phases in certain regions of phase space because of their high entropy of mixing,3 low curvature energy, and low packing frustration.2 The three commonly observed cubic surfaces in naturally occurring lyotropic amphiphilic systems are the primitive (P phase, Im3m space group), diamond (D phase, Pn3m space group), and gyroid (G phase, Ia3d space group) surfaces.2,4,5 These cubic phases find applications in materials science, in the synthesis of mesoporous nanomaterials with interesting structural and electronic properties, in biosensing, ultrafiltration, protein crystallization,4,6-8 and as chirally selective filters.5 Cellular membranes are composed of lipid molecules that can undergo self-assembly into a range of mesophases. Cubic phases are found in intracellular organelles such as the endoplasmic reticulum, Golgi apparatus, and mitochondria.2,8-12 The cubic * To whom correspondence [email protected]. † E-mail: [email protected].

should

be

addressed.

E-mail:

phases are efficient space partitioners which provide a large active surface area and allow access for incoming and outgoing molecules across interfaces. They are also thought to play an important role as possible intermediates in various cellular processes including endocytosis, exocytosis, membrane budding, cellular adhesion, and fusion and digestion processes in the stomach.2,6-11 The pathways and kinetics of self-assembly and transitions between noncubic and cubic as well as cubic and cubic phases can provide important insight into cellular processes involving cubic phases and have recently been an active area of experimental investigation.7,12-14 Experimental studies of the selfassembly of one-, two-, and three-dimensional periodic structures in amphiphilic systems have provided a wealth of information about the phase diagrams as the amphiphile concentration, temperature, and pressure are varied. In amphiphilic systems, the periodic structure is sensitive to addition of oil, salts, and other amphiphilic species. An understanding of the conditions of formation and stability of these periodic mesostructures, as well as the ability to control their structural properties, has implications for their exploitation as functional materials. For example, in the application of cubic structures in cellular membranes for molecular recognition and controlled release drug delivery,7,10,14-16 properties such as the lattice parameter and molecular specificity become important. Mesophases in amphiphilic systems extend over length scales at the mesoscale (10-100 nm) and are generated by microscopic interactions between the molecules, thus making the study of their self-assembly a truly multiscale problem. Continuum Navier-Stokes equation based methods do not take into account particle interactions and hence are insufficient to accurately follow the interfacial dynamics, while microscopic molecular dynamics methods are too computationally expensive to follow the self-assembly processes on the time scales and length scales involved. Our approach is to use the lattice-Boltzmann method based on the Shan-Chen model17,18 for the interaction between

10.1021/jp0731506 CCC: $40.75 © 2008 American Chemical Society Published on Web 02/21/2008

Self-Assembly of Mesophases oil and water particles extended to include orientable amphiphilic dipoles,19,20 which describes both the interactions between the explicit amphiphile particles and oil/water particles and the amphiphile-amphiphile particle interactions. Our amphiphilic lattice-Boltzmann model has been used previously to describe amphiphilic system self-assembly dynamics and rheology. For binary amphiphile and water systems, the self-assembly of the primitive phase has been simulated,21 while in ternary systems, the self-assembly and dynamics of the gyroid and sponge phases have been studied in considerable detail.5,22,23 The role of finite size effects was studied in previous simulations of the self-assembly of the gyroid mesophase,5,24 and systems of size 643 lattice sites were found to exhibit essentially the same phenomenology as larger 1283 and 2563 systems. However, the larger (2563) systems, although self-assembled into the gyroid phase at identical system parameters as those of the smaller systems, were riddled with defects and contained multiple gyroid domains.5,25 The rheology of the ternary gyroid mesophase has also been studied in great detail and found to exhibit non-Newtonian viscoelastic behavior, as expected from its liquid-crystalline morphology.24,26,27 In particular, the gyroid phase exhibits enhanced viscosity as compared to the sponge and lamellar phases. Since an increase in amphiphile concentration in the lamellar phase does not confer a higher viscosity, the rigidity of the gyroid phase is attributed to its morphology rather than its high concentration of amphiphiles.24 An understanding of the origins of morphology becomes even more important if/when such functional properties arise purely from the structural characteristics of the mesophase. Here, we report lattice-Boltzmann simulations in which a ternary mixture of oil, water, and amphiphile self-assembles to form three-dimensional periodic cubic phasessthe primitive and diamond phases in addition to the previously reported gyroid phase.5 With this work, we demonstrate that the self-assembly of three of the most commonly observed cubic structures in nature can be simulated within a single kinetic model with shortrange interactions between the quasi-molecular species, without imposing a macroscopic free-energy functional governing largescale system dynamics. We also report results from parameter sweeps in which various types of mesophases self-assemble under the competing influence of oil-water repulsion and the amphiphile-amphiphile interactions. The sequence of mesophases formed with increasing oil-water repulsion and increasing system pressure corresponds to the experimentally determined phase diagram, further validating the suitability of our model to simulate the self-assembly of cubic phases and its potential use in the study of the transition kinetics and pathways in cubic and noncubic mesophases. We also present results for the range of mesophases observed upon varying the relative concentration of the amphiphilic species. The majority of the experimental investigations in the literature have focused on lipid-water mixtures, which assemble into bilayers whose midsurface lies on a triply periodic minimal surface. It has been noted28 that when the thickness of the bilayer is small compared to the lattice constant, it is appropriate to model them with single (monolayer) structures of a ternary system, where the water regions on both sides of the bilayer are distinguished artificially by labeling them “water I” and “water II”. Useful insight into the self-assembly of such systems can hence be gained from our ability to accurately simulate the self-assembly of monolayer structures of ternary amphiphilic fluids.

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2951 2. The Lattice-Boltzmann Model We use a multicomponent version of the lattice-Boltzmann method in which the oil, water, and amphiphiles are modeled explicitly as three independent species. For convenience, we refer to oil particles as being “red” (r) and water particles as being “blue” (b). The amphiphile molecules are modeled as possessing two different fragments (a blue “head” and a red “tail”), each having an affinity for one of the two immiscible components. The amphiphile molecules are assumed to carry a dipole vector (or “director”). This degree of freedom is further coarse-grained by introducing a dipole vector field d(x,t), which represents the average orientation of all amphiphilic particles present at site x at time t. The orientation of d is allowed to vary continuously in time, but its magnitude is constant. We set the magnitude of the dipole vector for each amphiphile molecule, d0, as unity.29 Incorporation of dipole-carrying amphiphilic fluid particles results in two new types of interactions (between amphiphilic and nonamphiphilic particles and among amphiphiles themselves), which depend not only on the relative distance between particles but also on the dipolar orientations. Assuming that the dipole head and tail have equal and opposite color “charge”, e ) (1 (where +1 is red and -1 is blue), D is the dimension of the lattice and only nearestneighbor interactions considered; the additional forces are given by the following equations20,30,31

Fσ,s(x,t) ) -2ψσ(x,t)gσs

( )

d(x + ci∆t,t)‚ I ∑ i*0

Fs,σ(x,t) ) 2ψs(x,t)d(x,t)‚

∑σ gσs ∑ i*0

cici D ψs(x + ci∆t,t) 2 ci (1)

( ) I-

cici D ψσ(x + ci∆t,t) (2) 2 ci

and

Fs,s(x,t) ) -

4D c

2

gssψs(x,t)

[ ] I-

∑ i*0

{

d(x + ci∆t,t)d(x,t) :

cici D ci + d(x + ci∆t,t)d(x,t)‚ci + c2i

}

d(x,t)d(x + ci∆t,t)‚ci ψs(x + ci∆t,t) (3) Here, ci (i ) 0, 1, ..., 18) are the nineteen discrete velocity vectors on the D3Q25 regular lattice and ∆t is the latticeBoltzmann time step. Fσ,s is the force acting on nonamphiphilic particles σ (water and oil) due to amphiphile dipoles, Fs,σ is the force acting on amphiphilic particles due to all nonamphiphilic particles, and Fs,s is the force among amphiphilic particles. The coupling constants gσs and gss determine, respectively, the strength of interactions between water/oil particles and amphiphiles and between amphiphilic particles. The summation in eqs 1-3 is over all discrete velocities ci except c0, which is of zero magnitude. The interactions between nonamphiphilic species are modeled as self-consistently generated mean field body forces

2952 J. Phys. Chem. B, Vol. 112, No. 10, 2008

Fσ (x,t) ≡ -ψσ(x,t)

Saksena and Coveney

ψσj (x′,t)(x′ - x) ∑σj gσ σj ∑ x′

where ψσ(x,t) is the so-called effective mass, which can have a general form for modeling various types of fluids; we use ψσ ) (1 - e-nσ), nσ being the number density of the nonamphiphilic species σ.17 The amphiphilic lattice-Boltzmann model has been discussed in more detail elsewhere,19,20,30,31 and here, we simply state the fundamental equations governing the dynamics of the fluid. The model is characterized by the following set of coupled equations30

f σi (x + ci∆t,t + ∆t) - f σi (x,t) ) -∆t

(

)

f σi - f σ(eq) i τ

σ

∑σj ∑j Λσσijj f σjj + ∑j Λσsij f sj

+

(5)

f si (x + ci∆t,t + ∆t) - f si (x,t) ) -∆t

( ) f si - f s(eq) i

+

s

τ

∑σ ∑j Λsσij f σj + ∑j Λssij f sj

(6)

τRFR nR

(7)

j Λσσ ij ) ωi

[

]

c ‚c 1 σ σj σ σj σ σj i j (δ c R c ) + R ci ‚aσ∆t i j 2 4 cs cs

(8)

[

(10)

]

d h (x,t) - d(eq)(x,t) τd

(11)

Here, d is the average dipole vector per site, whose relaxation is governed by the BGK eq 11. , f s(eq) , and d(eq) are suitably chosen In eqs 5, 6, and 11, f σ(eq) i i local equilibrium distribution functions;31 τd is the relaxation time for the dipole vector, and d h is the average dipole at site x prior to collision. In our model for the amphiphilic species,30 the dipole vector distributions in eqs 10 and 11 are coupled to the effective mass ψs(x,t) via eqs 2, 3, 6, and 7. The dipole relaxation time τd in this model is taken to be a constant for simplicity. A more accurate but more computationally expensive description of amphiphilic dynamics would take into account the change in dipole relaxation time under the influence of the local fields. The diagonal components of the pressure tensor are computed as the sum of a kinetic term plus a virial mean field term24

∑R ∑k FRk (x)(ck - u(x))‚(ck - u(x)) + 1

gRRj ∑ [ψR(x)ψRj(x′) + ψRj(x)ψσ(x′)] × ∑ 4 R,Rj x′ (x - x′)‚(x - x′) (12)

The summation in the second term on the right-hand side runs over all combinations of the three fluid components, oil (r), water (b), and amphiphiles (s). Since the interaction matrix {gRRj} is symmetric with only nearest-neighbor interactions considered, the virial term reduces to

1 2

The τR is the relaxation time for the species R; this controls the kinematic viscosities of the corresponding fluid component 5,21-24 and is set to unity here for all fluid components. j σs sσ ss Additionally, the terms Λσσ ij , Λij , Λij , and Λij are matrix elements of the collision operators which result from mean field interactions among different fluid components.30 For example, j Λσσ ij is given by

∑i f si (x - ci∆t,t)d(x - ci∆t,t)

d(x,t + ∆t) - d h (x,t) ) -∆t

P(x) ≡

Equations 5 and 6 describe the time evolution of discrete velocity distribution functions f σi and f si belonging to component σ (oil and water) and amphiphile (s), respectively, and with discrete velocity ci. The first term on the right-hand sides of eqs 5 and 6 is the standard BGK collision operator.29,32 The effect of the interaction forces on the collisional relaxation is incorporated by shifting the average velocity used to calculate the local equilibrium distribution functions, f R(eq) , where index i R ) (σ, s) runs over all three distinct species present in the fluid, in the BGK collision operator by19

δuR )

f s(x,t)d h (x,t) )

(4)

gRRj ∑ ψR(x)ψRj(x + ck)ck‚ck ∑ R,R j k

(13)

Equation 12 is used to calculate the pressure of the fluid system. From eq 12, we can change the system pressure by changing the value of the interaction parameters, gss, gbr, and gbs. In particular, increasing the oil-water repulsion by making gbr more positive, keeping all other interaction parameters and species’ densities constant, causes the system to reach steady state at higher pressure, while a more negative value of the amphiphile-amphiphile coupling parameter, gss, keeping all other model parameters constant, causes the system to reach steady state at a lower pressure within the range of systems studied in this work.

with 3. Simulations σ j

F Rσ σj )

n

σ j

σ σ j

n

×

τ

Fσj

(9)

∑σj σj τ

Here, cs is the speed of sound which, in the D3Q25 lattice scheme we have used in this work, is equal to 1/x3, and aσ ) Fσ/Fσ is the force per unit mass. Hydrodynamic quantities such as the velocities uR of each fluid component are obtained from velocity moments of the corresponding distribution functions. The time evolution of amphiphilic dipoles d(x,t) is given by

We performed parameter sweeps by varying the (i) oil-water repulsion parameter gbr, (ii) amphiphilic coupling parameter gss, and (iii) relative concentration of oil, water, and amphiphilic species. The range of interaction parameters used in these simulations was limited by the numerical stability of the latticeBoltzmann algorithm and available computational resources. As described in the previous section, the simulation parameter gbr governs the interaction between the nonamphiphilic water (blue, σ ) b) and oil (red, σ ) r) species. A positive gbr is used to model repulsion between the species.19,20 The amphiphilic coupling parameter gss controls the interaction between amphiphile particles. For gss > 0, the coupling parameter models repulsion between two head groups (and attraction between a

Self-Assembly of Mesophases head and a tail), while gss < 0 describes attraction between two head groups (and repulsion between a head and a tail). A negative value of the water-amphiphile coupling parameter, gbs < 0, favors assembly of structures in which the head groups of the amphiphilic molecules are oriented toward the waterrich region and the tail groups are oriented toward the oil-rich region. We have used a constant value of gbs ) -0.006 throughout this work. The simulation cell is comprised of 643 lattice sites and is initialized in a homogeneous state in which the single particle distribution functions, fR(0) (where R can be oil (r), water (b), or amphiphile (s) particles), are distributed uniformly between the corresponding concentration parameter FR, that is, 0 e fR(0) e FR‚FR is specified with the constraint FR e 1.0. The system is then allowed to evolve under the collision-advection scheme of the lattice-Boltzmann method in the presence of interaction forces as described in the previous section eqs 5-11. Pressure is calculated as per eqs 12 and 13 and reported here in lattice units. The unit of mass in lattice units is the mass of a single particle; as in all of our previous work, we have taken the mass of oil, water, and amphiphile particles to be equal to unity. The unit of length is the lattice spacing and the unit of time is one lattice-Boltzmann time step. These lattice units have been used to report pressure in previous rheological simulations performed with this code;24,26,27 in particular, the pressure calculation has been used to predict the viscoelastic behavior of the gyroid mesophase, including relaxation times for the switch-over from solid-like to liquid-like behavior under oscillatory shear.27 Pressure values calculated from double-precision checkpoints are averaged over the 643 lattice sites and over 1000 time steps at steady state. The lattice-Boltzmann method provides a powerful approach for capturing the universal features of nonequilibrium fluid behavior. In the case of simple fluids, correspondence with experiment may be achieved by matching various dimensionless groups (e.g., Reynolds, Schmidt, Peclet numbers, etc.). For complex fluids such as those under investigation here, a complete matching to laboratory units cannot be easily achieved.33 The magnitude of the spatial and temporal lattice units can be estimated from comparisons with experimental data. Typical values for the lattice spacing and time step are 2 nm and 10 µs, respectively, based on observed space and time scales during self-assembly of amphiphilic cubic mesophases.22 Likewise, the mass can be equated with that of around six surfactant molecules, that is, 10-24 kg, based on observations of micelles formed from clusters of particles around a single lattice site, which approximate the aggregation number of monomers in a typical micelle containing ∼50 surfactant molecules. We find the widest range of mesophases with minimal change in interaction parameters for the parameter set Fr ) 0.3, Fb ) 0.3, and Fs ) 0.9, gbr ) -0.01, gbs ) -0.006, and variable gss. This parameter set was chosen to most efficiently demonstrate, with limited available computational resources, that the latticeBoltzmann method is able to exhibit self-assembly of the mesophases, in agreement with experiment in certain ranges of the parameter space. This parameter set thus forms the basis of Tables 1 and 2, where we additionally vary gbr to higher values, 0.05, 0.08, and 0.1, in order to access higher pressures of the system at steady state. In Table 3, the values of gss, gbr, and gbs have been chosen based on previous simulations using this LB3D code which resulted in the formation of the gyroid phase for Fr ) 0.7, Fb ) 0.7, and Fs ) 0.9.5 We vary the relative concentrations of the amphiphilic species for the state points in Table 3 to show

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2953 TABLE 1: Self-Assembly of Mesophases for Parameter Set 1: Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.0005, gbs ) -0.006, and Varying gbr, Displayed as a Function of Decreasing Pressure. The System Undergoes Coarsening for Increasing Values of the Repulsion Parameter gbr

mesophase description

(Pxx + Pyy + Pzz)/3

0.1 0.08

lamellar lamellar with periodic undulations lamellar with periodic undulations primitive

0.2571 0.2549

0.05 0.01

0.2517 0.2472

TABLE 2: Self-Assembly of Mesophases for Parameter Set 2: Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gbs ) -0.006, and Variable gbr and gss, Displayed as a Function of Decreasing Pressure gbr

gss

mesophase description

(Pxx + Pyy + Pzz)/3

0.1 0.05

-0.0005 -0.0005

0.2571 0.2517

0.01 0.01 0.01

-0.0005 -0.001 -0.002

0.01 0.01

-0.0025 -0.003

0.01

-0.011

lamellar lamellar with periodic undulations primitive hexagonal lamellar with periodic undulations lamellar lamellar with sponge defects sponge-like phase with lamellar domains

0.2472 0.2468 0.2458 0.2453 0.2447 0.2328

TABLE 3: Self-Assembly of Mesophases for Parameter Set 3 with Different Relative Concentrations of Oil, Water, and Amphiphile Species, gss ) -0.003, gbr ) 0.08, and gbs ) -0.006, Displayed as a Function of Decreasing Pressure Fr

Fb

Fs

mesophase description

0.9 0.7 0.5 0.4

0.9 0.7 0.5 0.4

0.9 0.9 0.9 0.9

sponge gyroid lamellar lamellar

(Pxx + Pyy + Pzz)/3 0.4715 0.4055 0.3304 0.2909

TABLE 4: Self-Assembly of Mesophases for Parameter Set 4 with Different Relative Concentrations of Oil, Water, and Amphiphile Species, gss ) -0.0005, gbr ) 0.01, and gbs ) -0.006, Displayed as a Function of Decreasing Pressure Fr

Fb

Fs

mesophase description

(Pxx + Pyy + Pzz)/3

0.9

0.9

0.9

0.4480

0.6

0.6

0.9

0.5

0.5

0.9

0.3

0.3

0.9

lamellar with periodic convolutions lamellar with periodic convolutions metastable state with hexagonal, primitive, and lamellar domains primitive

0.3469 0.3135 0.2472

the sequence of formation of sponge f gyroid f lamellar with increasing concentration of the amphiphilic species, in agreement with experiments and idealized phase diagrams.1,16 Finally, in Table 4, we explore the effect of variation in relative amphiphile concentration starting from the parameter set which gives the primitive phase in Table 2 and find a transition from the cubic (primitive) phase to the lamellar phase upon increasing the concentration of the nonamphiphilic species, in contrast to the trend reported in Table 3. The most interesting self-assembly behavior observed in previous work has been sponge and gyroid formation for a high amphiphile fraction of Fs ) 0.9 and is one of the key starting points in our parameter searches. Because of the inherent limitations of the lattice-Boltzmann scheme’s approach to equilibrium, there is no a priori reason to assume that for lower

2954 J. Phys. Chem. B, Vol. 112, No. 10, 2008 amphiphile concentrations, Fs, the species’ densities, Fb, Fr, and Fs, may not independently affect the self-assembly behavior, that is, act as three separate degrees of freedom instead of one compositional degree of freedom. Hence, in this work, we report interesting self-assembly behavior in a limited region of the model’s parameter space with Fs ) 0.9 and make comparisons with experiment. Historically,5,34 binary and ternary mixtures symmetric in oil and water densities have been simulated with this latticeBoltzmann code. This choice was initially made to study phase separation in binary mixtures, which was found to proceed via a spinodal decomposition path instead of the more complex nucleation-growth mechanism for the asymmetric mixture.34 We have adhered to ternary compositions which are symmetric in the oil and water densities in the present study to minimize complexity in the parameter space exploration. As has been mentioned in previous related work,34 the approach to equilibrium in most lattice-Boltzmann models is not guaranteed, owing to the lack of an H theorem, and in the case of interacting, multicomponent lattice-Boltzmann models, there is no well-established thermohydrodynamic theory stipulating equilibria to which the automaton can relax. Hence, lattice-BGK stationary regimes should be treated with caution and contrasted with experiment. It should be noted here that experimentally, self-assembly processes in amphiphilic mixtures often exhibit metastability and hysteresis,1,35 and we cannot rule out the possibility that intermediate metastable phases are a result of actual physical processes being modeled rather than a simulation artifact. The objective of this work is to identify regions of the model’s parameter space where we are able to see mesophase self-assembly and compare such behavior with experiment. As mentioned previously, our parameter searches have been limited by the stability of the algorithm and computational resources available. Previous investigations into lattice effects for the gyroid phase5 have confirmed that the orientation of the self-assembled gyroid is not influenced by the lattice, and we do not expect these to play a significant role in the current investigations. The alignment of the 643 lattice site systems in these simulations with the orthogonal coordinate axes of the simulation cell is a result of periodic boundary conditions which are used along the three coordinate directions. This effect is also observed in off-lattice dissipative particle dynamics simulations36,37 with periodic boundary conditions in which a lamellar monodomain aligned nontrivially with respect to the simulation cell is formed from an initial random mixture. Upon going to larger system sizes, we expect the influence of the periodic boundary conditions on the self-assembly to be mitigated, as in our previous work on the gyroid phase.5 Numerical instabilities are found to occur in our simulations of ternary amphiphilic fluids due to body force terms introduced in the lattice-Boltzmann equation to model interactions between the component species. These instabilities occur when the force coupling constants gbr, gbs, and gss or the mean densities of the fluid species are increased beyond certain threshold values. For these system parameters, the forcing term can become large enough to make the right-hand side of the lattice-Boltzmann equation (eqs 5 and 6) negative,22 which, in turn, leads to the values of the system properties becoming larger than can be represented as a double floating-point number on available machine architectures. In exploring the parameter space in the simulations reported here, we have reached such threshold values of the interaction parameters at different densities of the component species, where the lattice-Boltzmann equation becomes unstable. In addition to the obvious generation of

Saksena and Coveney

Figure 1. Simulated self-assembly of a ternary amphiphilic mixture into the gyroid phase on a 643 lattice; simulation snapshot taken after 15000 time steps.5 This is an isosurface representation showing the interface between the oil (red) and water (blue) domains.

spurious values for physical properties, we can detect the onset of instability when the particle velocities start approaching the speed of sound on the D3Q25 lattice, the value of which is 1/ x3. The lattice-Boltzmann simulations were performed on HPCx38 in the U.K. and across various multiprocessor machines on the U.S. production grid infrastructure, TeraGrid.39 We used the recently developed grid middleware known as the Application Hosting Environment (AHE)40 to launch simulations on HPCx and the TeraGrid machines. AHE simplifies the management of simulations, including the transfer of input and output files to grid resources. AHE provides a single, uniform interface to all of the grid resources and has built-in provenance functionality, allowing the user to store and review simulation workflow and histories, including pointers to simulation input and output. This makes the AHE a powerful tool for keeping track of simulations and exploring systems with complex parameter spaces. 4. Results We have identified the self-assembly of the primitive phase and the diamond phase in these lattice-Boltzmann simulations. The gyroid phase (see Figure 1) had been previously located using the same model for amphiphilic systems.5 Figure 1 is an isosurface representation of the self-assembled gyroid system showing the interface between the oil (red) and water (blue) domains. The interface between oil and water is defined as the surface on which the color order parameter, φ(r) ) nr(r) - nb(r) is equal to zero, where nr,b(r) is the number density of the red and blue particles, respectively, at point r in the system. This isosurface representation is used to elucidate the morphology of self-assembled structures reported here. The self-assembled primitive phase obtained from the evolution of an initially random ternary mixture under the lattice-Boltzmann collisionadvection scheme on a 643 lattice is shown in Figure 2a. For a clearer view, a 163 sublattice has been extracted from the simulation data set and is rendered in Figure 2b. The periodic nodal surface approximation of the primitive surface is given by FP ) cos(2πx) + cos(2πy) + cos(2πz).7 An artificially generated 163 sublattice, which was extracted from the periodic nodal surface approximation corresponding to the simulated primitive surface, is shown in Figure 2c for comparison. Figure 2(d) contains a common depiction of the primitive surface.5,7

Self-Assembly of Mesophases

Figure 2. (a) Simulated self-assembly of a 643 lattice site system into the primitive phase; simulation snapshot taken after 15000 time steps. Simulation parameters are Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.0005, gbr ) 0.01, and gbs ) -0.006. (b) A 163 sublattice extracted from the simulated system in Figure 2. (c) 163 sublattice extracted from an artificially generated periodic nodal representation of the primitive surface with 19 unit cells on a 643 lattice; see text. (d) An artificially generated primitive surface with three unit cells on a 643 lattice.

The surface in Figure 2d is a periodic nodal representation containing three unit cells in a 643 lattice system. In comparison, the simulated primitive system in Figure 2a contains approximately 19 unit cells per 643 lattice system, resulting in a lower-resolution primitive surface with a smaller lattice parameter. The self-assembled diamond phase and the artificially constructed surface from the periodic nodal surface approximation, FD ) cos(2πx)cos(2πy)cos(2πz) - sin (2πx)sin(2πy)sin(2πz),7 are illustrated in Figure 3a and b, respectively. A magnified 163 sublattice from the simulated diamond phase is shown in Figure 3c, with the corresponding periodic nodal surface approximation in Figure 3d. We are also able to locate the two-dimensional periodic hexagonal phase (H); the 643 simulated lattice system selfassembled into the hexagonal phase and a 163 sublattice extracted from the simulated system are illustrated in Figure 4a and b, respectively. We have also located a one-dimensional periodic lamellar phase (L) and another lamellar phase with periodic convolution, which are shown in Figures 5 and 6, respectively. Next, we report simulations in which self-assembly of a range of mesophases is observed as the system parameters are systematically varied. In parameter set 1 (see Table 1), the simulation parameters Fr, Fb, Fs, gss, and gbs are kept constant, and the oil-water repulsion parameter, gbr, is varied. Increasing gbr corresponds to the formation of mesophases in the order primitive f lamellar with periodic undulations f lamellar. As the repulsion between the oil and water species increases with increasing gbr, the system is expected to be driven toward an energetically favorable domain-segregated arrangement with the amphiphiles largely located at the interface between domains of substantially increased characteristic size. We observe such

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2955

Figure 3. (a) Simulated self-assembly of a 643 lattice site system into the diamond phase; simulation snapshot taken after 15000 time steps. Simulation parameters are Fr ) 0.05, Fb ) 0.05, Fs ) 0.7, gss ) 0.00035, gbr ) 0.01, and gbs ) -0.006. (b) Artificially generated diamond surface with 10 unit cells. (c) A 163 sublattice extracted from the simulated system in Figure 3, which is self-assembled into the diamond phase. (d) A 163 sublattice extracted from the artificially generated diamond surface in Figure 3b.

Figure 4. (a) Simulated self-assembly of a 643 lattice site system into a hexagonal phase; simulation snapshot taken after 15000 time steps. Simulation parameters Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.001, gbr ) 0.01, and gbs ) -0.006. (b) A 163 sublattice extracted from the simulated system in Figure 4.

a coarsening of domains with increasing gbr, which is also accompanied by an increase in the pressure of the fluid system (eq 12). In all three lamellar systems in Table 1, the 643 system self-assembles into alternating 16 lamellae of oil and water and thus have the same interlamellar spacing. However, the intensity of the peak in the structure factor versus wavenumber profile, which is a measure of the characteristic domain mass, increases with gbr, thus confirming that the system becomes increasingly segregated with increasing repulsion. The peak in the plot for all three values of gbr occurs at wavenumber k ) 1.6 (inverse lattice units), and this corresponds to the repeat period of the water and oil density modulations. The structure factor is defined here as the Fourier transform of the color order parameter φ(r).5,34

2956 J. Phys. Chem. B, Vol. 112, No. 10, 2008

Figure 5. Simulated self-assembly of a 643 lattice site system into a lamellar phase; simulation snapshot taken after 15000 time steps. Simulation parameters are Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.0005, gbr ) 0.1, and gbs ) -0.006.

Figure 6. Simulated self-assembly of a 643 lattice site system into a lamellar phase with periodic convolutions on the lamellar surface; simulation snapshot after 15000 time steps. Simulation parameters are Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.0005, gbr ) 0.05, and gbs ) -0.006.

Using parameter set 2 (see Table 2), by varying the colorcolor repulsion parameter gbr and amphiphile-amphiphile coupling parameter gss, we are able to induce the self-assembly of a range of mesophases. The order of appearance of mesophases with decreasing system pressure is lamellar f lamellar with periodic convolutions f primitive f hexagonal f lamellar f sponge-like phase. The sponge-like phase is shown in Figure 7; it is a random microemulsion and contains interpenetrating regions of oil and water. The order of appearance of the mesophases for parameter sets 1 and 2 is in accordance with the pressure jump experiments of Seddon et al.13 and Czeslik et al.35 in which the mesophases are observed in the sequence lamellar f primitive f hexagonal with decreasing pressure. However, in our simulations, we also observe a second lamellar phase and a coarsened sponge-like phase upon going to even lower pressures. For parameter set 3 (see Table 3), with gbr ) 0.08 and gss ) -0.003, upon lowering the oil + water concentration (Fb + Fr) or, equivalently, increasing the relative amphiphile concentration (Fs), the system self-assembles into a lamellar state. Upon decreasing the relative amphiphile concentration (Fs), the system

Saksena and Coveney

Figure 7. Simulated self-assembly of a 643 lattice site system into a sponge-like phase with coarsened domains; simulation snapshot taken after 15000 time steps. Simulation parameters are Fr ) 0.3, Fb ) 0.3, Fs ) 0.9, gss ) -0.011, gbr ) 0.01, and gbs ) -0.006.

Figure 8. Simulated self-assembly of a 643 lattice site system into a bicontinuous sponge phase; simulation snapshot taken at 15000 time steps. Simulation parameters are Fr ) 0.7, Fb ) 0.7, Fs ) 0.4, gss ) -0.003, gbr ) 0.08, and gbs ) -0.006.

self-assembles into a gyroid phase; further decreasing of the amphiphile concentration results in a bicontinuous sponge phase (see Figure 8). This corresponds to the experimentally observed phase sequence for cubic lipid bilayer systems15 and previously published idealized phase diagrams,1,16 where decreasing amphiphile concentration corresponds to a transition from a lamellar phase to a cubic mesophase. Interestingly, for parameter set 4 (see Table 4), gbr ) 0.01 and gss ) -0.0005, we see the reverse trend, in which increasing the relative amphiphile concentration favors the cubic primitive phase. Such deviations from ideal behavior have been observed in experiments for certain lipid-water mixtures; cubic phases have also been observed on the water-poor side of the lamellar phases for binary aqueous solutions of the anionic soap sodium diethylhexylsulfosuccinate and binary aqueous solutions of galactolipids and lecithins.1 For the state point Fr ) 0.5, Fb ) 0.5, and Fs ) 0.9 in Table 4, a metastable state is observed. The stable state for this state point is marginal between the hexagonal and lamellar phase, with periodic convolutions depending on the starting configuration of the system.

Self-Assembly of Mesophases 5. Conclusion The interactions between the various components of amphiphilic mixtures play an important role in determining the mesophases obtained as a result of self-assembly. In this paper, we report results on the influence of the oil-water repulsion parameter, amphiphile-amphiphile coupling parameter, and relative amphiphile concentration on self-assembly in ternary mixtures simulated using the lattice-Boltzmann method. We are able to simulate a wide range of mesomorphic behavior in agreement with experimental studies of lyotropic solutions. In certain regions of parameter space, we also observe deviations from ideal phase behavior3,16 as the concentration of the components is varied. Experimentally, such deviations have been observed for amphiphilic systems.16 Some of our future work will focus on determining the extent of finite size effects and the role of hysteresis in amphiphilic self-assembly processes. Hysteresis, whereby a metastable intermediate phase can persist on the order of days, weeks, or months, has been observed in experiments.1,35 In our latticeBoltzmann simulations, however, it is possible to perturb the system out of a metastable state by changing the quasi-molecular interaction parameters or by imposing nonequilibrium boundary conditions such as external shear. Performing parameter sweeps for larger systems, for example, comprised of 1283 and 2563 lattice sites, in order to determine phase behavior as well as the role of finite size effects and hysteresis is extremely computationally expensive; we plan to pursue these issues in the future. The role of high-performance computational grids is important in these studies in order to harness the power of multiple highperformance resources for the rapid execution of computational steering-type parameter sweeps5,41 and geographically distributed simulations,42,43 which are very effective in reducing the time needed to obtain results. In conclusion, we have been able to simulate the primitive, diamond, and gyroid cubic mesophases as well as the noncubic hexagonal and lamellar mesophases using our bottom-up latticeBoltzmann method. The ability to simulate such rich phase behavior with the lattice-Boltzmann method provides us with the necessary framework within which to conduct further equilibrium and nonequilibrium studies of noncubic T cubic phase and cubic T cubic phase-transition dynamics. It is also the starting point for a systematic study of the rheological properties of these mesophases. Acknowledgment. This research has been funded by EPSRC through the RealityGrid (GR/R67699), RealityGrid Platform (EP/C536452/1), Large Scale Lattice-Boltzmann Simulation of Liquid Crystals (EP/E045111/1), and Rapid Prototyping of Usable Grid Middleware (GR/T27488/01) Grants and through the OMII Robust Application Hosting using WSRF::Lite (GR/ 290843/01) Project. Access to the U.S. TeraGrid was provided undertheNRACandPACSGrantsMCA04N014andASC030006P. Data transfer from supercomputing resources in the U.K. and U.S. was facilitated by UKLight within the EPSRC ESLEA Grant GR/T04465. References and Notes (1) Fontell, K. Colloid Polym. Sci. 1990, 268, 264. (2) Winter, R.; Erbes, J.; Czeslik, C.; Gabke, A. J. Phys.: Condens. Matter 1998, 10, 11499.

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2957 (3) Safran, S. A. In Statistical Thermodynamics of Surfaces, Interfaces and Membranes; Pines, D., Ed.; Frontiers in Physics. Perseus Books Group: New York, 1994. (4) Schwarz, U. S.; Gompper, G. In Morphology of Condensed Matter: Physics and Geometry of Spatially Complex Systems; Mecke, K. R., Stoyan, D., Eds.; Springer: New York, 2002; Vol. 600 of Lecture Notes in Physics, pp 107-151. (5) Chin, J.; Coveney, P. V. Proc. R. Soc. London, Ser. A 2006, 462, 3575; 10.1098/rspa.2006.1741. (6) Squires, A. M.; Templer, R. H.; Seddon, J. M.; Woenckhaus, J.; Winter, R.; Finet, S.; Theyencheri, N. Langmuir 2002, 18, 7384. (7) Squires, A. M.; Templer, R. H.; Seddon, J. M.; Woenkhaus, J.; Winter, R.; Narayanan, T.; Finet, S. Phys. ReV. E 2005, 72, 011502. (8) Marrink, S. J.; Tieleman, D. P. Biophys. J 2002, 83, 2386. (9) Schwarz, U. S.; Gompper, G. Langmuir 2001, 17, 2084. (10) Marrink, S. J.; Tieleman, D. P. J. Am. Chem. Soc. 2001, 123, 12383. (11) Schwarz, U. S.; Gompper, G. J. Chem. Phys. 2000, 112, 3792. (12) Conn, C. E.; Ces, O.; Mulet, X.; Finet, S.; Winter, R.; Seddon, J. M.; Templer, R. H. Phys. ReV. Lett. 2006, 96, 108102. (13) Seddon, J. M.; Squires, A. M.; Conn, C. E.; Ces, O.; Heron, A. J.; Mulet, X.; Shearman, G. C.; Templer, R. H. Philos. Trans. R. Soc. London, Ser. A 2006, 364, 2635. (14) Angius, R.; Murgia, S.; Berti, D.; Baglioni, P.; Monduzzi, M. J. Phys.: Condens. Matter 2006, 18, S2203. (15) Shearman, G. C.; Ces, O.; Templer, R. H.; Seddon, J. M. J. Phys.: Condens. Matter 2006, 18, S1105. (16) Kaasgaard, T.; Drummond, C. J. Phys. Chem. Chem. Phys. 2006, 8, 4957. (17) Shan, X.; Chen. H. Phys. ReV. E 1993, 47, 1815. (18) Shan, X.; Chen, H. Phys. ReV. E 1994, 49, 2941. (19) Chen, H.; Boghosian, B. M.; Coveney, P. V.; Nekovee, M. Proc. R. Soc. London, Ser. A 2000, 456, 2043. (20) Love, P. J.; Nekovee, M.; Coveney, P. V.; Chin, J.; GonzalezSegredo, N.; Martin, J. M. R. Comput. Phys. Commum. 2003, 153, 340. (21) Nekovee, M.; Coveney, P. V. J. Am. Chem. Soc. 2001, 123, 12380. (22) Gonza´lez-Segredo, N.; Coveney, P. V. Europhys. Lett. 2004, 65, 795. (23) Gonza´lez-Segredo, N.; Coveney, P. V. Phys. ReV. E 2004, 69, 061501. (24) Gonzalez-Segredo, N.; Harting, J.; Giupponi, G.; Coveney, P. V. Phys. ReV. E 2006, 73, 031503/1. (25) Harting, J.; Harvey, M. J.; Chin, J.; Coveney, P. V. Comput. Phys. Commun. 2004, 165, 97. (26) Giupponi, G.; Coveney, P. V. Math. Comput. Simul. 2006, 72, 124. (27) Giupponi, G.; Harting, J.; Coveney, P. V. Europhys. Lett. 2006, 73, 533. (28) Schwarz, U. S.; Gompper, G. Phys. ReV. E 1999, 59, 5528. (29) Chen, S.; Doolen, G. Annu. ReV. Fluid Mech. 1998, 30, 329. (30) Nekovee, M.; Coveney, P. V.; Chen, H.; Boghosian, B. M. Phys. ReV. E 2000, 62, 8282. (31) Chen, H.; Boghosian, B. M.; Coveney, P. V.; Nekovee, M. Proc. R. Soc. London, Ser. A 2000, 456, 2043. (32) Bhatnagar, P. L.; Gross, E. P.; Krook, M. Phys. ReV. 1954, 94, 511. (33) Cates, M. E.; Desplat, J. C.; Stansell, P.; Wagner, A. J.; Stratford, K.; Adhikari, R.; Pagonabarraga, I. Philos. Trans. R. Soc. London, Ser. A 2005, 363, 1917. (34) Gonzalez-Segredo, N.; Nekovee, M.; Coveney, P. V. Phys. ReV. E 2003, 67, 046304(17). (35) Czeslik, C.; Winter, R.; Rapp, G.; Bartels, K. Biophys. J. 1995, 68, 1423. (36) Guo, H.; Kremer, K.; Soddemann, T. Phys. ReV. E 2002, 66, 061503/1. (37) Guo, H.; Kremer, K. J. Chem. Phys. 2003, 119, 9308. (38) HPCx Home Page. http://www.hpcx.ac.uk. (39) TeraGrid. http://www.teragrid.org. (40) Coveney, P. V.; Saksena, R. S.; Zasada, S. J.; McKeown, M.; Pickles, S. Comput. Phys. Commun. 2007, 176, 406; 10.1016/ j.cpc.2006.11.011. (41) Chin, J.; Coveney, P. V.; Harting, J. Proceedings of the 3rd U.K. e-Science All Hands Meeting, Nottingham, U.K., Aug 31-Sept 2, 2004; Cox, S. J., Ed.; http://www.allhands.org.uk/2004/proceedings/papers/181.pdf. (42) Boghosian, B.; Finn, L. I.; Coveney, P. V. http://www.realitygrid.org/publications/GD3.pdf. (43) Karonis, N.; Toonen, B.; Foster, I. J. Parallel Distrib. Comput. 2003, 63, 551.