Theory and computer simulations of amphiphilic Janus particles · 2017-10-26 · We apply molecular...

Theory and computer simulations of

amphiphilic Janus particles

vorgelegt von

Diplom-Physiker

Gerald Rosenthal

aus Berlin

von der Fakultät II - Mathematik und Naturwissenschaften

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften

Dr. rer. nat.

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender: Prof. Dr. Martin Schoen

1. Gutachterin: Prof. Dr. Sabine H. L. Klapp

2. Gutachter: Dr. habil. Thomas Weikl

Tag der wissenschaftlichen Aussprache: 13.06.2012

Berlin 2012

D 83

Eidesstattliche Versicherung

Hiermit erkläre ich an Eides statt, dass ich die von mir vorgelegte Dissertation selbstständigangefertigt und die verwendeten Quellen und Hilfsmittel vollständig angegeben habe. DieKooperation mit anderen Wissenschaftlern habe ich in der Dissertation kenntlich gemacht.

Gerald Rosenthal

Abstract

This work is concerned with amphiphilic Janus particles. Nowadays Janus particles areexperimentally synthesizable and are of great importance in modern technological processes.Here, Janus particles are modeled as (hard or soft) spheres composed of two distinct surfaces,one hydrophilic and the other hydrophobic. Their orientation is described by a vectorrepresenting an internal degree of freedom.

We apply molecular dynamics simulations and density functional theory to investigatethe structure formation of those amphiphilic Janus particles in the bulk phase, as well asin the vicinity of a planar wall or in confinement, such as a slit pore geometry. Our densityfunctional theory study of the bulk system focuses on determining an upper temperaturelimit for bilayer formation, where the reduced temperature is an inverse measure for theanisotropic coupling strength. Additionally, we discuss the possibility of a condensationphase transition.

A more detailed study of the bulk system is provided by analyzing energy fluctuationsand cluster size distributions in the framework of molecular dynamics simulations. Here wedetermine the aggregation line in a temperature-density-diagram. Below this aggregationline clusters with a size-dependency on density and reduced temperature are found. Theseclusters are mainly spherical, where we find a predominant role of icosahedrons in densersystems. Interestingly, we find no hint of a condensation transition of these clusteredsystems in our molecular dynamics simulations. A brief study on the incorporation ofhigher order terms in the pair interaction potential shows additional structures, such asbilayers and worm-like aggregates.

We additionally compute in our computer simulations dynamical quantities and corre-lation functions to achieve an insight in the time evolution of the system. Analyzing thetranslational mean-square displacement, we also observe indications of hindered diffusiondue to aggregation, whereby the bond autocorrelation function shows long-lived micellaraggregates.

Going back to the “basic” model we apply our density functional theory to study the com-petition between surface-induced structure formation and fluid-fluid interactions. Studyingdensity and polarization profiles as a function of the distance to the wall, we find surface-induced bilayer-like structure formation. Moreover, such structures can strongly be influ-enced by the hydrophobicity of the surface potential. The investigation of confinement ina slit-pore shows interesting frustration effects, which cause an oscillation of the normalpressure as a function of the wall distance. To compare the bulk structures with the struc-tures in confinement we additionally apply molecular dynamics simulations in the slit poregeometry. These reveal a stabilization of bilayers close to the walls for the dense system atsuitable temperatures, which indeed show a qualitative agreement with our density func-tional theory study. Our fundamental study of an elementary model of amphiphilic Janusparticles provides an interesting contribution to the general understanding of such Janusparticles and their aggregation behavior.

Zusammenfassung

Diese Arbeit widmet sich amphiphilen Janusteilchen. Janusteilchen sind heutzutage exper-imentell herstellbar und besitzen eine große Bedeutung in technologischen Prozessen. Indieser Arbeit werden die Janusteilchen als (harte oder weiche) Kugeln mit unterschiedlichenOberflächen modelliert, wobei die eine Halbkugel hydrophil und die andere entsprechend hy-drophob ist. Die Orientierung eines jeden Janusteilchens wird mittels eines Einheitsvektorsbeschrieben, der als innerer Freiheitsgrad fungiert.

Wir untersuchen die Strukturbildung dieser Teilchen mittels Molekulardynamiksimula-tionen (MD) und Dichtefunktionaltheorie (DFT) im unbegrenzten System, sowie in derGegenwart einer planaren Wand und in der räumlich begrenzten Geometrie einer Schlitz-pore. Mit Hilfe der DFT bestimmen wir eine obere Schranke für die Temperatur, unterwelcher sich Membranen bzw. Doppelschichten bilden. Die reduzierte Temperatur ist dabeiein inverses Maß für die anisotrope Kopplungsstärke. Darüberhinausgehend diskutieren wireinen möglichen Kondensationsphasenübergang.

Eine detailliertere Untersuchung des unbegrenzten Systems auf Basis von Energiefluktu-ationen und Clustergrößenverteilungen bietet sich im Rahmen der Molekulardynamiksim-ulation an. Unter Verwendung dieser bestimmen wir die Aggregationstemperaturlinie alsFunktion der Dichte. Ein Unterschreiten dieser Linie resultiert in der Formation von Clus-tern verschiedener Größe und ist abhängig von Dichte und Temperatur. Diese Cluster sindhauptsächlich sphärisch, wobei Ikosaeder bei hoher Dichte das System dominieren. DieMolekulardynamiksimulationen zeigen keinen Hinweis auf einen Kondensationsphasenüber-gang. Wir untersuchen ferner die Auswirkung eines Einbeziehens höherer Ordnungen in dasPaarwechselwirkungspotential, welches es ermöglicht Membranen und wurmartige Anord-nungen als weitere Strukturen zu erzeugen.

Neben der statischen Betrachtung des Systems berechnen wir in unseren Simulationendynamische Größen und Korrelationsfunktionen um über das zeitliche Verhalten Aussagentätigen zu können. Die Analyse der mittleren quadratischen Versetzung (mean-square-displacement) weist eine durch Aggregatbildung verzögerte Diffusion auf, wobei wir aus derAutokorrelation der Bindungen auf langlebige, mizellenartige Aggregate schließen.

Auf Grundlage des ursprünglichen Modells ohne höhere Ordnungen studieren wir überdiesmittels DFT den Wettbewerb von oberflächeninduzierter Strukturformation und derjenigendie im unbegrenzten System vorherrscht. Hierbei analysieren wir die Dichte- und Polarisa-tionsprofile mit Wanddistanzabhängigkeit, welche auf eine oberflächeninduzierte, membra-nartige Strukturbildung schließen lassen. Solche Strukturen können darüber hinaus mittelsder Hydrophobizität des Oberflächenpotentials kontrolliert und beeinflusst werden. Inter-essante wandabstandsabhängige Frustrationseffekte zeigen sich durch die beschränkendeGeometrie einer Schlitzpore, welche sich auch in einer Oszillation der Druckkomponentesenkrecht zu den Wänden ausdrückt. Zusätzliche MD-Simulationen von Janusteilchen ineiner Schlitzpore ermöglichen uns den Vergleich der Strukturbildung mit derjenigen im un-begrenzten System. Es zeigt sich bei großen Dichten und passenden Temperaturen eineStabilisierung von doppelschichtartigen Membranen in der Nähe der Wände. Dieses Ergeb-nis deckt sich qualitativ mit unserer Dichtefunktionaltheoriestudie. Unsere grundlegendeUntersuchung eines elementaren Modells für amphiphile Janusteilchen bietet einen Beitragzum generellen Verständnis dieser Janusteilchen und deren Aggregatbildung.

Publication list

• G. Rosenthal and S. H. L. Klapp, Ordering of Janus particles at planar walls: Adensity functional study, J. Chem. Phys. 134, 154707 (2011).

• G. Rosenthal, K. E. Gubbins, and S. H. L. Klapp, Self-assembly of model amphiphilicJanus particles, J. Chem. Phys. 136, 174901 (2012).

Contents

1. Introduction 1

2. Statistical mechanics background and the Hamiltonian 7

2.1. The Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2. Ensembles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.1. Microcanonical ensemble . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.2. Canonical ensemble . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2.3. Grand canonical ensemble . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3. Averages and the ergodic hypothesis . . . . . . . . . . . . . . . . . . . . . . 102.4. The factorization of the probability density . . . . . . . . . . . . . . . . . . 11

3. The underlying Model 133.1. The present model of amphiphilic Janus particles . . . . . . . . . . . . . . . 133.2. Model modifications for the molecular dynamics simulations . . . . . . . . . 173.3. Models for walls and confinement . . . . . . . . . . . . . . . . . . . . . . . . 18

4. Density functional theory 194.1. General introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194.2. Ideal gas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214.3. Fundamental measure theory . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4.3.1. Starting point of a fundamental measure theory . . . . . . . . . . . . 22

4.3.2. Scaled particle theory . . . . . . . . . . . . . . . . . . . . . . . . . . 254.3.3. The original Rosenfeld functional . . . . . . . . . . . . . . . . . . . . 274.3.4. White-Bear version mark II . . . . . . . . . . . . . . . . . . . . . . . 28

4.4. DFT approximations to handle pair interactions beyond a reference potential 294.4.1. Mean-field approximation . . . . . . . . . . . . . . . . . . . . . . . . 304.4.2. Modified mean-field approximation . . . . . . . . . . . . . . . . . . . 31

4.5. Confinement and walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.6. Density functional theory applied to our model . . . . . . . . . . . . . . . . 32

4.6.1. The homogeneous isotropic system in modified mean-field approxi-mation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.6.2. Parameters for our calculations . . . . . . . . . . . . . . . . . . . . . 34

5. Molecular dynamics simulations 355.1. The Verlet algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.2. Thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365.3. Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.4. Truncation and shifting of potentials . . . . . . . . . . . . . . . . . . . . . . 39

5.5. Characterizing equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.6. Finite size effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.7. Additional numerical details . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

6. Correlation functions and further relevant quantities 456.1. Structure correlation functions . . . . . . . . . . . . . . . . . . . . . . . . . 45

6.2. Cluster related order parameters and distributions . . . . . . . . . . . . . . 476.3. Time-dependent correlation functions . . . . . . . . . . . . . . . . . . . . . 486.4. Mean-square displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 496.5. Specific heat capacity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

6.6. Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

7. The bulk system of amphiphilic Janus particles 537.1. Density functional study of the bulk . . . . . . . . . . . . . . . . . . . . . . 537.2. Molecular dynamics study of the bulk . . . . . . . . . . . . . . . . . . . . . 56

7.2.1. Aggregation diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . 567.2.2. Cluster properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 597.2.3. Microscopic structure . . . . . . . . . . . . . . . . . . . . . . . . . . 617.2.4. Dynamic properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 657.2.5. The influence of the interaction range . . . . . . . . . . . . . . . . . 68

7.2.6. Results for the bulk system of modified Janus particles . . . . . . . . 69

8. The influences of surfaces on amphiphilic Janus particles 758.1. Density functional study of the influence of surfaces . . . . . . . . . . . . . 75

8.1.1. The neutral wall case . . . . . . . . . . . . . . . . . . . . . . . . . . 75

8.1.2. Influence of surface fields . . . . . . . . . . . . . . . . . . . . . . . . 788.1.3. Confinement between two neutral walls . . . . . . . . . . . . . . . . 81

8.2. Molecular dynamics study in confinement . . . . . . . . . . . . . . . . . . . 838.2.1. Microscopic structure . . . . . . . . . . . . . . . . . . . . . . . . . . 87

9. Conclusions and outlook 93

A. Appendix 99A.1. Yukawa potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99A.2. Soft wall potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

A.3. Maxwell construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101A.4. Reduced units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101A.5. Remark to the expansions of the pair potential and the pair correlation

function in rotational invariants . . . . . . . . . . . . . . . . . . . . . . . . . 102

B. Appendix 103

B.1. Proof for general introduction to DFT . . . . . . . . . . . . . . . . . . . . . 103B.1.1. Proof (1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103B.1.2. Proof (2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103B.1.3. Proof (3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

B.2. Numerical implementation of density functional theory . . . . . . . . . . . . 104

B.3. Functional derivation - how to do it . . . . . . . . . . . . . . . . . . . . . . 106B.4. Density functional theory calculations for Janus particles . . . . . . . . . . . 107

B.4.1. The anisotropic contribution to the grand canonical functional . . . 107B.4.2. Functional minimization with respect to α . . . . . . . . . . . . . . . 107B.4.3. Functional minimization with respect to ρ . . . . . . . . . . . . . . . 109

C. Appendix 113C.1. Cluster search algorithm and clustering criteria . . . . . . . . . . . . . . . . 113C.2. Forces and torques for molecular dynamics simulations of Janus particles . . 114

C.2.1. Basic Janus interaction . . . . . . . . . . . . . . . . . . . . . . . . . 114C.2.2. Incorporation of higher-order pair interactions . . . . . . . . . . . . 114C.2.3. Forces at an unstructured planar soft wall . . . . . . . . . . . . . . . 115

C.3. Cell lists and neighbor lists . . . . . . . . . . . . . . . . . . . . . . . . . . . 115C.4. Initial conditions for the molecular dynamics simulations . . . . . . . . . . . 115C.5. Conservation of total momentum in molecular dynamics simulations . . . . 115

D. Programs and libraries 117

1Chapter 1.

Introduction1

In our work we focus on the so-called Janus particles, which are subject of many recentstudies and applied in modern technologies, as we will see in the subsequent paragraphs.Here, the term Janus refers mythological to the two-faced Roman god, and nowadays spans awide class of particles composed of at least two chemically or physically distinctive surfaces.

One way to view Janus particles is that they represent a coarse-grained model of am-phiphilic molecules. The latter are generally composed of a polar, hydrophilic (i.e., “water-loving”) head and at least one hydrophobic (i.e., “fat-loving”) tail, typically a carbon chain.This coarse-graining treats such molecules as (condensed to) a sphere with anisotropic sur-face properties. The presence of these two ingredients, in combination with an aqueoussolvent, yields a large variety of self-assembled structures on different length scales. Theserange from molecular size micelles to mesoscopic membranes, bicontinuous foams, andlamellar phases. From an application point of view, amphiphiles are used in a variety ofcontexts, e.g., to reduce the surface tension in complex mixtures (such as water and oil) likeemulsifier in food, tertiary crude-oil recovery [3], or as detergents. The impact of surfacesis also important in a variety of applications, including synthesis of nanoporous materials[4], thin film deposition for lithographic processes and devices [5], enhancement of chemicalreactions [6], and, more recently, to stabilize bundle and network formation in solutions ofcarbon nanotubes [7, 8].

From the theoretical side, the self-assembly of amphiphilic molecules has been investi-gated by a variety of approaches and models [9, 10, 11], including lattice gas systems [12],Ginzburg-Landau theory [10], density functional studies of entropic models (where the sur-factant is represented by a sphere plus an infinitely thin rod [13, 14, 15], or by a dimer [16]),and off-lattice simulations [17, 18] of flexible bead-spring molecules [19], with and withoutexplicit solvent.

But Janus particles are more than just a coarse-grained model of surfactants or moregeneral amphiphilic molecules. In 1989 Casagrande et al. [20] introduced the term Janusparticles to describe coated glass spheres with amphiphilic properties in the sense that one ofthe hemispheres was hydrophilic (attracting water) and the other one hydrophobic (avoidingwater). The hydrophilicity was maintained on one side by a protecting cellulose varnish,when the other side was prepared with octadecyltrichlorosilane (OTS) which enables acoverage of aliphatic chains, which give the hydrophobic properties. Since then, and asanticipated by de Gennes on the occasion of his Nobel lecture [21], the surface propertiesof these Janus particles have become an area of great interest on their own. Irrespective ofthe chemical details, the directional dependence of the interactions between Janus particlestypically leads to rich phase behavior including self-assembly into aggregates such as micelles

1Selections from this chapter have been reprinted with permission from our studies [1, 2]. Copyright 2011and 2012, respectively, American Institute of Physics.

1

1. Introduction

Figure 1.1.: Some techniques of synthesizing Janus particles sketched similar to Ref. [32], a)masking technique, b) use of reactive fluxes or fields, c) microcontact printing,d) contact with reactive medium.

[22, 23, 24, 25, 26] or chains [27], as well as formation of mesoscale structures (such asmembranes) [25, 26, 28].

Significant experimental progress has been made in synthesizing such particles, e.g., bymasking techniques, techniques using a reactive directional flux or field, microcontact print-ing, and other processes involving a reactive medium [22, 29, 30, 31, 32, 33] since then. Wesketch these four techniques in Fig. 1.1. This variety of manufacturing processes has alsosubstantially widened the range of possible properties of Janus particles and the character ofthe resulting (anisotropic) interactions between such particles. For example, several recentexperimental studies have suggested Janus particles with electric dipolar or quadrupolar[27, 34], or magnetic dipolar character [22, 32, 35, 36].

A wide range of recent applications and technological processes incorporate Janus par-ticles. Electronic paper [37, 38] and sensors [39] use electronic Janus particles. Further,magnetic Janus particles can be used as microrheological probes [40, 41]. Of course the“original” amphiphilic Janus particles are still important, as they can, e.g., serve as stabi-lizers in emulsions [42, 43] and modify textiles such that they become water-repellent [44].Moreover, there is a strong fundamental interest triggered by both experiments and theoryto understand these self-assembly processes of Janus-like particles as a bottom-up processfor the design of future nanomaterials [45].

For bulk systems scanning electron microscopy [23] is the method of choice to experimen-tally capture detailed information on self-assembled structures. Self-assembled structuresin the vicinity of surfaces can be experimentally studied using, e.g., neutron reflectometry(targeting the thickness of the layer) [46], grazing incidence small-angle neutron scatter-ing [47], and atomic force microscopy (targeting the lateral structure) [48], or by studying

2

Figure 1.2.: Sketch of a Janus particle, a) one-patch particle [25], b) present work andoriginally from Ref. [28]. The core region of the Janus particle is inpenetrablein both models. Here, the left side (orange) is hydrophobic (attractive) in a)and b), where in b) the right side (blue) is additionally hydrophilic (repulsive).Contrary to a) are the interaction regions in b) not homogeneous, i.e., theinteraction strength is strongest at the left and right sides and decreasing withthe distance from the particle as well as to the top and bottom.

adsorption isotherms [49].

One of the first models to characterize and study amphiphilic Janus particles from atheoretical point of view has been proposed by Tarazona and coworkers [28] in 1995. Tara-zona and coworkers originally followed the idea of a coarse-grained model of amphiphilicmolecules. Their model treats the solvent implicitly and describes the amphiphilic moleculesas spherical particles composed of two hemispheres, one hydrophilic and the other hydropho-bic. The solvent-mediated interaction is taken into account through an effective directionalpotential involving orientation vectors defined by the symmetry axis of the spheres. Thispotential results from low order contributions gained from expanding a general pair inter-action potential of linear molecules in rotational invariants. Despite these simplifications(which completely neglect geometrical factors, such as relative size of the head group), themodel is capable of describing bilayer, vesicle, and micelle formation. We sketch such aparticle in Fig. 1.2b).

There are other independently derived models [22, 50, 51] for Janus particles, which arein fact quite similar to the model introduced by Tarazona and coworkers [28] in their generaldescription. Probably the most popular model goes back to the one-patch model suggestedby Kern and Frenkel [52] and has recently been studied by Sciortino et al. [25]. This“patchy-particle” model involves localized surface areas exerting homogeneous attractiveforces, whereby the core is impenetrable. In Fig. 1.2 we show a comparison of the Janusparticles studied by Sciortino et al. a) and those considered by Tarazona and coworkersb). In fact, the “patchy-particle” model was subject of extensive Monte Carlo computersimulations and integral equation theory studies [25, 26, 53], in which micelles, vesicles andbilayers were found. In addition to those various aggregation structures, the simulationsreveal the existence of a condensation transition, below which the system separates in a

3

1. Introduction

Figure 1.3.: The figure shows experimental epifluorescence microscopy results (green back-ground) and Monte Carlo simulation results of Janus particles, where the yellowside is hydrophobic and the blue side is negatively charged. The different pic-tures show different salt concentrations: A) deionized Janus particles in water,B) 1mM KNO3, C) 5mM KNO3. The Janus particles with a diameter of1µm are created by coating the later hydrophobic hemisphere with gold andfunctionalization of the same side with octadecanethiol after the coating. Thehydrophilic side is negatively charged due to carboxylic acid groups. Reprinted(adapted) from Ref. [22] with permission [Copyright (2012) American ChemicalSociety].

dilute phase of micelles, and a denser phase of larger clusters. The corresponding phasecoexistence curve becomes narrower and shifted toward higher densities upon cooling. Themodel studied with Monte Carlo simulations in Ref. [22] describes the Janus particles ashydrophobic on one side and charged on the other. Granick and coworkers [22] use a veryshort range interaction for the attractive hydrophobic sides, while they apply three differentrepulsive interactions with a Yukawa-like distance dependency resulting from the negativelycharged hydrophilic side. Granick and coworkers compare their simulation results withthose of Luijten and coworkers [22], we show in Fig. 1.3 a reprint of their results, whichshow micelles [cf. Fig. 1.3 A), B) and D)] and worm-like structures [cf. Fig. 1.3 C)]. Butstill a general understanding of the role of the various terms (and their prefactors) in therotationally invariant expansion on the overall self-assembly and phase behavior is missing.

Contrary to the extensively discussed case of true surfactant molecules, the self-assemblyof Janus particles at surfaces has so far been rarely considered, an exception being a recentstudy of Hirose et al. [54], who used a macroscopic theory based on the Young’s equation.

Against this background, the aim of this thesis is twofold. First, we study amphiphilicJanus particles in the bulk system. Later, we incorporate planar surfaces and surfacepotentials preferring specific orientations of the Janus particles. We also consider (strongly)confined geometry like slit pores. The Janus particles are described following the originalwork of Tarazona [28], where we mainly consider the “basic” model involving only the lowestorder terms of a general pair interaction of linear molecules.

Our choice of this particular model is stimulated by the density functional theory calcu-

4

lations in Ref. [28], which show the capability of the “basic” model to describe micelle aswell as bilayer formation. Interestingly, micelles are indeed the most important structuresof real amphiphilic Janus particles in the experimentally relevant density range[22, 23].

In our work we apply two different methods. We start with a classical equilibrium densityfunctional theory (DFT) study [55, 56], where we mainly employ a mean-field approxima-tion to treat the anisotropic part of the interactions, while the repulsive [hard sphere (HS)]contribution is treated on a more sophisticated level. Specifically, we employ the so-calledfundamental measure theory (FMT) (Ref. [57]) which has turned out to be extremely suc-cessful for the description of inhomogeneous hard sphere systems [58, 59, 60]. Densityfunctional theory provides a relatively fast method to deal with large systems (comparedto molecular dynamics and Monte Carlo simulations) and is originally formulated in thegrand canonical ensemble. Moreover, inhomogeneous systems such as the influence of sur-faces can be studied. Within this method our investigations focus on planar structures,such as bilayers. The key question is to determine thermodynamic conditions under whichself-assembly arises at the surface, in contrast to the corresponding bulk system. More-over, we explore the impact of different surface properties concerning their hydrophilic orhydrophobic character in particular. Second, we apply equilibrium molecular dynamics(MD) simulations to study the aggregation behavior and associated dynamics of the Janusmodel, which is related to considerations of the bulk within the framework of DFT. Atthis point, we find some advantages of MD simulations with respect to DFT. Indeed, onemain technical restriction of the DFT approach (beyond the unavoidable approximation fortreating the interparticle correlations) is that it requires an input for the shape of expectedaggregates, such as a (perfectly spherical) micelle or a (perfectly planar) bilayer. Clearly,this restriction is absent in a computer simulation. The comparison to a Monte Carlomethod shows the advantage of studying easily time evolution in MD simulations. On theother hand, we had to introduce some small model modifications, which go along with aMD study, such as a soft core repulsion instead of the hard core. Further, we extended theanisotropic potential to distances of possible overlap of two particles. Here we study themodel for a broad range of temperatures and densities within the fluid phase regime. Ourresults show that below a density-dependent aggregation temperature the (model) Janusparticles self-assemble into a rather wide variety of cluster types. Moreover, particularly atintermediate densities the aggregation process strongly affects the translational and rota-tional dynamics. Lastly, we extend the MD simulations to study confined geometries andcompare these results with those given through our DFT treatment.

This thesis is organized as follows. We start with a brief introduction of the generalHamiltonian and mention some relevant statistical ensembles in Ch. 2. The underlyingmodel of amphiphilic Janus particles is presented in Ch. 3. Our methods are described inthe subsequent chapters, where we derive the density functional theory in Ch. 4 and themolecular dynamics simulation in Ch. 5. Further, we introduce correlation functions, orderparameters and other important properties in Ch. 6. We split our results in two parts, wherewe focus on the bulk system in Ch. 7 and on the competition between surface-fluid and fluid-fluid interactions in Ch. 8. The conclusions and outlook follow in Ch. 9. Moreover, we addan elaborate appendix, where we show more details of various derivations and numericalimplementation.

5

1. Introduction

6

2Chapter 2.

Statistical mechanics background andthe Hamiltonian

In this chapter we briefly introduce the underlying Hamiltonian in Sec. 2.1 and considerthe important statistical ensembles of our work in Sec. 2.2. Further, we briefly commenton averages in Sec. 2.3, which are of general importance for both the derivation of thedensity functional theory and our molecular dynamics simulation results. Finally we facethe factorization of the probability density to find the points of the phase space in Sec. 2.4.

2.1. The Hamiltonian

A good starting point for our theoretical considerations is the Hamiltonian H, which givesthe energy of a system depending on the positions, momenta, orientations and angularmomenta of the particles in the system. The time dependency of the trajectories is notconsidered as long as ensemble averages are taken into account, as done in our classicalequilibrium density functional theory. On the other hand the trajectories in our moleculardynamics simulations depend on the time. We follow arguments of Gubbins [60] in thissection and assume that the Hamiltonian splits into a classical and a quantum mechanicalpart. A first approximation is that both parts do not interact, which is a good approxi-mation for equilibrium quantities. The quantum mechanical part describes vibrations andinternal rotations. In our study we concentrate on coarse-grained particles and thereforeomit the quantum mechanical contribution. We are then left with the classical part whichis given as

H(

rN ,pN ,pNω ,ω

N)

=N∑

i=1

p2i

2mi+

N∑

i=1

pTωi

· Ii · pωi

2+ Epot

(

rN ,ωN)

, (2.1)

where ri is the position of the particle i, pi is its momentum, ωi = (φi, θi, ψi)T is a vec-

tor giving the three Euler angles, φi, θi, ψi, and pωiis the conjugate momentum to ωi.

The potential energy Epot

(

rN ,ωN)

splits into the contribution of the particle interactions

U(

rN ,ωN)

and the contribution of the external potentials V(

rN ,ωN)

. The mass is given

by mi, Ii is the tensor that describes the moment of inertia and T gives the transposition

of a vector. The contribution U(

rN ,ωN)

simplifies for pair interactions, φ (ri, rj ,ωi,ωj),to

U(

rN ,ωN)

=1

2

N∑

i=1

N∑

j=1,j 6=i

φ (ri, rj ,ωi,ωj) . (2.2)

7

2. Statistical mechanics background and the Hamiltonian

In the case of linear molecules ωi reduces to (ϕ, θ) with ϕ = 0 . . . 2π and θ = 0 . . . π. If Ii

can be diagonalized we can choose a rotation to obtain the body-fixed principal axes α andthen write the rotational kinetic contribution as a function of Jiα, the components of theangular momentum with respect to the axis α. Using the relation Ji = Iipωi

we find

N∑

i=1

pTωi

· Ii · pωi

2=

N∑

i=1

JTi · I−1

i · Ji

2=

N∑

i=1

JTi R · diag

(

I−1i

)

RT · Ji

2=

N∑

i=1

∑

α=x,y,z

J2iα

2Iα,

(2.3)

where the rotated angular momentum Jiα depends on both ωi and pωi, R is a rotation

matrix and Iα is the eigenvalue corresponding to axis α of the diagonalized Ii. At this pointwe insert Eq. (2.3) in Eq. (2.1) and obtain

H(

rN ,pN ,pNω ,ω

N)

=N∑

i=1

p2i

2mi+

N∑

i=1

∑

α=x,y,z

J2iα

2Iα+ Epot

(

rN ,ωN)

. (2.4)

We note, that in the molecular dynamics simulations the angular velocities pωiare renamed

as ωi.

2.2. Ensembles

2.2.1. Microcanonical ensemble

The microcanonical ensemble [61] describes a system which is isolated. In fact, moleculardynamics simulations are performed in this ensemble, if no thermostat or other constraintsare applied. The number of particles N and the volume V are fixed within a system inthe microcanonical ensemble. Further, the energy of a realistic system can contrary to amodel system vary slightly between E and E + ∆, where ∆ is small with respect to E.Additionally the total momentum and the total angular momentum are zero. As a startingpoint we consider the probability density fmc to find the points of the phase space withinthe energy shell E . . . E + ∆, that is

fmc

(

rN ,pN ,pNω ,ω

N)

=

1

Γ (E) ∆, E ≤ H

(

rN ,pN ,pNω ,ω

N)

≤ E + ∆

0, otherwise, (2.5)

where Γ (E) ∆ is the volume of the energy shell. This equation gives a constant probabilityfor every point of the phase space within the energy shell. We now consider an infinite thinenergy shell

fmc

(

rN ,pN ,pNω ,ω

N)

≡ lim∆→0

fmc

(

rN ,pN ,pNω ,ω

N)

(2.6)

=Γ (E)−1 δ(

E −H(

rN ,pN ,pNω ,ω

N))

, (2.7)

where δ (. . .) is the Dirac δ-distribution. The probability density is normalized as follows

Trmcf(

rN ,pN ,pNω ,ω

N)

≡ 1. (2.8)

8

2.2. Ensembles

Taking this normalization into account and inserting Eq. (2.7) we can derive

Γ (E) =1

hnNN !

∫∫∫∫

d3Nrd3Npd3Npωd3Nω′δ

(

E −(

H(

rN ,pN ,pNω ,ω

N)))

(2.9)

=Trmcδ(

E −(

H(

rN ,pN ,pNω ,ω

N)))

. (2.10)

In the equation above is h the Planck constant and n gives the degrees of freedom of aparticle. The integrations are carried out with respect to the positions

d3Nr = ΠNi=1d

3ri = ΠNi=1Πj=x,y,zdrij , (2.11)

the momenta

d3Np = ΠNi=1d

3pi = ΠNi=1Πj=x,y,zdpij , (2.12)

the angular velocities

d3Npω = ΠNi=1d

3pωi= ΠN

i=1Πj=x,y,zdpωij(2.13)

and the Euler angles

d3Nω′ = ΠNi=1dφidθidψi (2.14)

of each particle (molecule) i.

2.2.2. Canonical ensemble

We introduce the canonical ensemble with analogous considerations and definitions to themicrocanonical ensemble. The canonical ensemble [61] describes a system which is in con-trast to the microcanonical ensemble in contact with a heat bath, which fixes the temper-ature T , i.e. energy fluctuations occur. Therefore, the canonical ensemble is described byconstant V , N , T . An important quantity for averages is the probability density,

f(

rN ,pN ,pNω ,ω

N)

= ZN (T, V )−1 e−β(H(rN ,pN ,pNω

,ωN)), (2.15)

where β = (kBT )−1 with the temperature T and the Boltzmann constant kB. The normal-ization condition for the probability density,

Trcf(

rN ,pN ,pNω ,ω

N)

≡ 1, (2.16)

leads us to the partition function,

ZN (T, V ) =1

hnNN !

∫∫∫∫

d3Nrd3Npd3Npωd3Nω′e−β(H(rN ,pN ,pN

ω,ωN)) (2.17)

=Trce−β(H(rN ,pN ,pN

ω,ωN)). (2.18)

Equation (2.18) provides us with an expression for the Helmholtz free energy,

F (T, V,N) = −β−1lnZN (T, V ) . (2.19)

As mentioned in Sec. 2.2.1 molecular dynamics simulations are native in the microcanonicalensemble. It is tempting to control the temperature in molecular dynamics simulations bya thermostat to achieve a simulation in the canonical ensemble, where we can fix thetemperature, but it is not an easy task and some thermostats do not reach this goal orpossess other limitations (see Sec. 5.2).

9


2.2.3. Grand canonical ensemble

The grand canonical ensemble [61] allows energy and particle number fluctuations and isdescribed by fixed T , V and chemical potential µ. This ensemble will be our choice for thederivation of the density functional theory (DFT) in Sec. 4. In the DFT we will need theprobability density in the grand canonical ensemble, that is

fN

(

rN ,pN ,pNω ,ω

N)

= ΞN (T, V )−1 e−β(HN(rN ,pN ,pNω

,ωN)−µN), (2.20)

with the normalization condition

TrgcfN

(

rN ,pN ,pNω ,ω

N)

≡ 1. (2.21)

The grand partition function Ξµ (T, V ) occurring in Eq. (2.20) follows from the normaliza-tion condition as

Ξµ (T, V ) =∞∑

N=0

1

hnNN !

∫∫∫∫

d3Nrd3Npd3Npωd3Nω′e−β(HN(rN ,pN ,pN

ω,ωN)−µN)

(2.22)

=Trgce−β(HN(rN ,pN ,pN

ω,ωN)−µN). (2.23)

With these quantities we obtain an expression for the grand canonical potential, whichreads

Ω (T, V, µ) = −β−1lnΞµ (T, V ) = −PV, (2.24)

where P is the pressure.

2.3. Averages and the ergodic hypothesis

A statistical description of a system in equilibrium holds only for averages, where we havetwo different ways to define those. First, we consider the ensemble average of a quantity Qin thermodynamic equilibrium, which is defined as follows

〈Q〉 = Tr[

f(

rN ,pN ,pNω ,ω

N)

Q(

rN ,pN ,pNω ,ω

N)]

, (2.25)

where N gives the number of particles. Such an average will be important for our DFTstudy. The second one is the time average, which is defined as follows

〈Q〉t =1

t

∫ t

0Q(

rN ,pN ,pNω ,ω

N)

dt′, (2.26)

where for finite time t the result depends on the initial conditions. Therefore, the existenceand independence from such initial conditions is postulated for limt→∞ 〈Q〉t. Both, theensemble and the time averages are assumed to be equal in the ergodic hypothesis, i.e. whenthe trajectory of a particle goes “through” every point of the phase space and is independentof the initial conditions.

10

2.4. The factorization of the probability density

2.4. The factorization of the probability density

Our aim is to carry out the integration of the momenta and the angular momenta in statis-tical quantities, like the partition function. Our steps in this section follow the derivationin Ref. [60]. We know from Eqs. (2.3) and (2.4), that the rotational kinetic contribution tothe Hamiltonian depends on JN , which is in general determined by both pN

ω and ωN . Thepotential energy depends also on ωN . Therefore, we have to consider the following generalinequality for the probability density [cf. Eqs. (2.7), (2.15), and (2.7)]

f(

rN ,pN ,pNω ,ω

N)

6= f(

pN)

f(

pNω

)

f(

rN ,ωN)

. (2.27)

To overcome this problem, we will switch our description from the variable set(

pN ,ωN)

to(

JN ,ωN)

. By considering the normalization condition

∫∫∫∫

d3Nrd3Npd3Npωd3Nω′f

(

rN ,pN ,pNω ,ω

N)

= 1, (2.28)

where d3Nω′ is defined by Eq. (2.14), we can define an analogous normalization conditionwith the new variable set:

∫∫∫∫

d3Nrd3Npd3NJd3Nω′f ′(

rN ,pN ,JN ,ωN)

= 1. (2.29)

Here we use the definition d3NJ = ΠNi=1d

3Ji = ΠNi=1Πj=x,y,zdJij and we perform the trans-

formation of the variable set [60]

f ′(

rN ,pN ,JN ,ωN)

= f(

rN ,pN ,pNω ,ω

N)

N∏

i=1

sinθi. (2.30)

At this point we define

f(

rN ,pN ,JN ,ωN)

=f ′(

rN ,pN ,JN ,ωN)

∏Ni=1 sinθi

(2.31)

and

d3Nω = d3Nω′N∏

i=1

sinθi. (2.32)

With these two definitions and the transformation to the new variable set we can conclude

f(

rN ,pN ,JN ,ωN)

= f(

pN)

f(

JN)

f(

rN ,ωN)

. (2.33)

From now on, the trace Tr will refer to the new variable set, where J is used instead of pω

and where we will integrate over d3Nω instead of the primed expression.

11


12

3Chapter 3.

The underlying Model

This chapter introduces the present model of amphiphilic Janus particles. The model itselfand its origin are found in Sec. 3.1 and some minor modifications for our molecular dynamicssimulations are emphasized in Sec. 3.2. This chapter closes with our treatment of walls andconfinement in Sec. 3.3.

3.1. The present model of amphiphilic Janus particles2

In this study, we employ a simple, coarse-grained model of an amphiphilic system originallysuggested by Tarazona and coworkers [28]. The solvent, which is omnipresent in real am-phiphilic solutions, is treated implicitly. The derivation starts from linear molecules. Forthose it is possible to expand the pair interaction potential in rotational invariants [28, 60]in a space fixed frame (coordinate system), yielding

φ (r,ω1,ω2) =∑

l1,l2,l

∑

m1,m2,m

ξl1l2,l (r)C (l1, l2, l;m1,m2,m)

· Yl1,m1 (ω1)Yl2,m2 (ω2)Y ∗l,m (ω) , (3.1)

where C (l1, l2, l;m1,m2,m) are Clebsch-Gordon coefficients and Yl,m (ω) is a spherical har-monic, with ∗ giving the complex conjugate. Further, ξl1l2,l (r) are radial-dependent ex-pansion coefficients. Due to rotational symmetry we can rotate the system such that theintermolecular axis r = r1 − r2 = r12 becomes parallel to the z-axis by the choosingω = (0, φ), as shown in Ref. [60]. This rotation leads to an intermolecular frame with therotated molecule orientations ω′

1 and ω′2. The change into the intermolecular frame yields

Y ∗l,m (0, ϕ) =

√

2l + 1

4πδm0, (3.2)

with δm0 = 1, if m = 0 and δm0 = 0 otherwise.

At this point we go into more specific details. Since we are interested in Janus particlesrather than molecules or surfactants in the present model, the amphiphilic Janus particlesare represented by spherical particles consisting of two hemispheres, one hydrophilic andthe other hydrophobic. We introduce u1 and u2, which are unit vectors denoting theorientations of the spheres. Specifically, ui points from the hydrophobic to the hydrophilicside of particle i. The vectors ui = u (ω′

i) (i = 1, 2) are parameterized by ω′i = (ϕi, θi),

with ϕi = [0, 2π] and θi = [0, π]. With these unit vectors, the unit vector rji = −rij/rij

2Selections from this section have been reprinted with permission from our studies [1, 2]. Copyright 2011and 2012, respectively, American Institute of Physics.

13

3. The underlying Model

Figure 3.1.: Total pair potential for three different configurations.

and rij = |rij| we can replace cosθ1 = u1 · r21 and cosθ2 = −u2 · r21. At this point we canwrite Eq. (3.1) in the intermolecular frame as follows

φ(

r12,ω′1,ω

′2

)

=∑

l1,l2,m

ξl1l2,m (r)Yl1,m

(

u1(

ω′1

)

· r21, ϕ1)

Yl2,−m

(

−u2(

ω′2

)

· r21, ϕ2)

,

(3.3)

with

ξl1l2,m (r) =∑

l

√

2l + 1

4πC (l1, l2, l;m,−m, 0) ξl1l2,l (r) . (3.4)

The corresponding low-order spherical harmonics can be found in the Appx. A.5. In thelowest order we get an isotropic contribution, which will be associated with the hard-sphere(HS) repulsion ξ0,0,0 (r12) ∝ φHS (r12) to create an impenetrable particle core. An attractiveisotropic interaction would cause a segregation of solvent and Janus particles, which wouldnot make sense in our understanding of amphiphilic Janus particles. In the next order wedefine ξ1,0,0 (r12) = ξ0,1,0 (r12) ∝ φ1 (r12), because of symmetry reasons, i.e. interchange-ability of the particles. The proportionality relation is used to omit constants arising fromthe spherical harmonics in the later derivations. With this low-order approximation of ageneral pair interaction potential of linear molecules we see, that the total pair potentialbetween two Janus-like particles at positions r1 and r2 subdivides into a hard-sphere (HS)and an anisotropic contribution,

φ (r12, u1, u2) = φHS (r12) + φI (r12, u1, u2) , (3.5)

where the HS potential is given by

φHS (r12) =

∞, r12 < σ,

0, r12 ≥ σ, (3.6)

14

3.1. The present model of amphiphilic Janus particles

with σ being the HS diameter. The anisotropic interaction is defined as

φI (r12, u1, u2) = φ1 (r12) (u1 − u2) · r21, (3.7)

where

φ1 (r12) =

0, r12 < σ,

C exp (−λ (r12 − σ)) /r12 r12 ≥ σ(3.8)

is a Yukawa potential, which is, e.g., used to describe the screened interactions of electriccharges, see Appx. A.1. In Eq. (3.8) the parameters C and λ measure the coupling strengthand the (inverse) range of the interaction, respectively. If not stated otherwise, we setλ∗ = λσ−3 = 3.

In Fig. 3.1 the radial dependency of the pair potential [cf. Eq. (3.5)] for some relevantconfigurations of two particles is shown. Within our model the particles prefer to be ori-entated in opposite directions, such that the hydrophobic sides point toward one another.The opposite configuration (with facing hydrophilic sides) is energetically least favorable.Parallel orientations are energetically neutral. The strong energetic preference (or penalty)of configurations with facing hydrophobic (hydrophilic) sides mimics the effects expectedin a real system which would include, e.g., water as a solvent. Clearly, the water moleculeswill preferentially adsorb at the hydrophilic side of each particle. The resulting steric ex-clusion yields an effective repulsion of the hydrophilic sides of neighboring Janus particles.On the other hand, the fact that the hydrophobic sides dislike water effectively favors con-figurations where these sides point toward each other, such that the contact with water isminimized. An additional source of attraction between hydrophobic sides can arise from afunctionalization with, e.g., thiol. An example has been described in Ref. [23] where thiol(and other molecular groups) is used to functionalize the gold-coated side of a polystyrenesphere. The functionalization yields a depletion area with respect to water and thus, anadditional effective attraction between the hydrophobic (gold) hemispheres.

Moreover, the basic model for our amphiphilic Janus particles can be extended by thenext order of Eq. (3.3), which yields

φI (r12, u1, u2) =φ1 (r12) (u1 − u2) · r21 + q1φ2 (r12)Q1 (r21, u1, u2)

+ q2φ3 (r12)Q2 (r21, u1, u2) , (3.9)

where

Q1 (r21, u1, u2) = − (u1 · r21) (u2 · r21)

−√

1 − (u1 · r21)2√

1 − (u2 · r21)2cos (ϕ1 − ϕ2)

= − u1 · u2, (3.10)

and

Q2 (r21, u1, u2) =1

2

(

3 (u1 · r21)2 + 3 (u2 · r21)2 − 2)

. (3.11)

15


For simplicity we have defined ξ1,1,0 (r12) ∝ q1φ1 (r12) and ξ2,0,0 (r12) = ξ0,2,0 (r12) ∝q2φ1 (r12). In fact, Tarazona and coworkers [28] also included the last term of Eq. (3.9)in their pair interaction, where the aim was to stabilize parallel orientations, such as mem-branes. The middle term was omitted in there study, because its incorporation wouldenforce a numerical procedure to calculate α (r,ω) in Eq. (4.15), which otherwise can bedone analytically. We limit our density functional theory study to the pair interactionpotential given in Eq. (3.5), which includes the most important properties of amphiphilicJanus particles. In our molecular dynamics simulations we study both the pair interactionpotential defined by Eq. (3.5) and the pair potential incorporating additional interactionsas given in Eq. (3.9) with q1 = 0 and q2 = 0.5. The results for the latter potential arediscussed in Sec. 7.2.6. We note, that some small modifications of the pair interactionpotential are made in our MD simulation study, which are discussed in Sec. 3.2.

It is worthwhile and interesting to consider three of the models mentioned in the intro-duction and to briefly compare the present model with those. The potential suggested byHess and coworkers [50] contains, in addition to the angle-dependent terms appearing al-ready in Eq. (3.5), terms involving Eq. (3.10), as well as higher powers of the various scalarproducts. The general form of these terms can also be rationalized on the basis of thegeneral expansion of the interaction potential between two linear molecules into sphericalharmonics [50, 60] as done in our case. However, contrary to our model favoring antiparallelalignment (see Fig. 3.1), the prefactors of the various terms in Ref. [50] were chosen suchthat parallel orientation of the Janus spheres is preferred.

Another Janus potential has been suggested by Sciortino et al. [25]. In this model, the(spherical) particles have one homogeneous patch; the patches of neighboring particles theninteract via a square-well potential. The remaining parts of the spheres only induce stericrepulsion. This picture is quite different to our model, which should rather be comparedto a particle with two patches, one mimicking the hydrophilic part and the other one thehydrophobic part. Specifically, in our model the strongest attractive (repulsive) interactionsoccur when the particles are aligned in an antiparallel manner with facing hydrophobic(hydrophilic) sides. Contrary to that, the one-patch model in Ref. [25] concentrates onthe attraction arising between the particles if the two patches become coupled through thesquare-well zone.

The model used by Granick and coworkers [22] applies contrary to us different interactionstrength and interaction range for the hydrophilic and hydrophobic sides of the Janusparticles. The hydrophobic hemispheres attract one another on very short distances oftwo particle centers not exceeding 1.05σ. This attraction is not distance dependent butit changes with the relative position of both particles. The hydrophilic side is modeled ascharged, where a Yukawa potential is used to describe the distance dependency. The authorsdiffer between three different repulsions, where each is homogeneous. These three differentcases are: Hydrophilic hemispheres facing toward one another, head-to-tail configurationsand hydrophobic sides facing toward one another, the latter at a distance of the particlecenters of at least 1.05σ. The major difference to our model is, that we describe thehydrophilic and hydrophobic side by the same model. Moreover, we do not introduce arepulsion of head-to-tail configurations as long as we do not consider extensions of ourmodel as done in Eq. (3.9).

16

3.2. Model modifications for the molecular dynamics simulations

Figure 3.2.: Total pair potential for different configurations, kBT/ǫ = 1 and kBTσ/C =0.167. Reprinted with permission from our study [2]. Copyright 2012, Ameri-can Institute of Physics.

3.2. Model modifications for the molecular dynamics simulations

Hard spheres are not easy to implement in a standard molecular dynamics simulation.In fact, particle collision events must be traced and then, when such an event occurs themomenta have to be redirected as done in an event driven simulation. Therefore, we use softspheres instead of hard spheres in our molecular dynamics simulations. Our pair interactionpotential is then

φ (r12, u1, u2) = φsoft (r12) + φI (r12, u1, u2) , (3.12)

where the soft-sphere potential is given by

φsoft (r12) = 4ǫ (σ/r12)12 . (3.13)

Further we redefine Eq. (3.8) appearing in Eqs. (3.7), (3.9) and (3.12) as follows

φ1 (r12) = Cexp (−λ (r12 − σ)) /r12, (3.14)

because the hard core is not present anymore and particles can interpenetrate each other.In fact, the molecular dynamics simulations show a penetration depth into a sphere notexceeding ≈ 6% of σ within our parameters. The radial dependency of Eq. (3.12) for threedifferent configurations is plotted in Fig. 3.2.

In our simulations, we need an explicit moment of inertia for the Janus particles. Herewe take the one of a solid sphere, so

I =

∫ 2π

0dϕ

∫ π

0dΘsinΘ

∫ R

0drρr2(rsinΘ)2 =

2mR2

5. (3.15)

17


3.3. Models for walls and confinement

Models of confinement are manifold, reaching from structureless walls to the incorporationof the inner composition. One of the simplest models is a structureless, planar hard wallwith the surface potential

φsurf (z) =

∞, z < σ/2,

0, z ≥ σ/2. (3.16)

We use such a wall potential in our density functional study. To include effects of preferentialadsorption in our DFT study we also consider potentials of the form

φsurf (z, u) =

∞, z < σ/2,

φsurf0 (z) · u, z ≥ σ/2

, (3.17)

where we use the expression

φsurf0 (z) = C2 exp (−λ2 (z − σ/2)) ez/z. (3.18)

It is also useful to define

φsurf (z) = φsurf0 (z) Θ

(

z − σ

2

)

(3.19)

for our later derivations. In the Eq. (3.18) the (unit) vector ez points along the z-directionand u [cf. Eq. (3.17)] is the orientation of the Janus particle considered. The sign ofthe parameter C2 then determines which side of the particle is attracted by the surface.Specifically, negative (positive) values of C2 correspond to a preference of the hydrophobic(hydrophilic) side. Finally, the range of the surface potential is controlled by the parameterλ2. For simplicity, we set λ∗

2 = λ2σ3 = λ∗. In experiments silicon surfaces can be coated

with polymer films of varying thickness to vary the hydrophobicity [62]. Of course the lackof an explicit solvent in our model causes an attraction between our hydrophilic surfaceand the Janus particles, where in an experiment the solvent would cover the surface andthe Janus particles would probably not really be in contact with the surface.

In the molecular dynamics simulations we describe a wall according to

φsoft wall (x) =4

45ρwallπǫwall

σ12

|x− xwall|9, (3.20)

which is an unstructured, planar soft-wall as derived in Appx. A.2. We chose a dense wall,characterized by ρwallσ

3ǫwall = ǫ.

18

4Chapter 4.

Density functional theory

Density functional theory (DFT) [63, 64, 65] is a widely used method in physics and chem-istry. It is used for quantum mechanical systems, like systems considered in the solid statephysics, as well as in the soft-matter physics. As mentioned in Sec. 2.1 quantum mechanicalinfluences can be neglected in our study. Therefore, it is useful to consider the classicalDFT equivalent, which will be used in this work.

We start with the general introduction of DFT in Sec. 4.1 and then consider the simplestcase of the ideal gas in Sec. 4.2. The next more complicated system, which is often usedas a reference systems, is composed of hard spheres. We recapitulate the derivation of themodern incorporation of hard spheres, the so-called fundamental measure theory, in Sec. 4.3.To include more general pair interactions, we take a look at the so-called λ-expansion inSec. 4.4. After that, in Sec. 4.5, some considerations of confinement and walls are made andwe close this chapter with the DFT applied to our model of amphiphilic Janus particles inSec. 4.6.

4.1. General introduction

Our general introduction follows closely Evans [55]. First, we consider the grand canonicalensemble, which is the native ensemble for density functional theory - but it is worthmentioning that representations in the canonical ensemble also exist [66]. The equilibriumprobability density f0 [see Eq. (2.20)] is our starting point:

f0 = Ξ−1e−β(HN −µN) = fN , (4.1)

where HN [see Eq. (2.1)] describes the Hamiltonian in a system of N particles, Ξ is thegrand partition function, and µ denotes the chemical potential. We now consider thegeneral functional

Ω [f ] = Tr[

f(

HN − µN + β−1lnf)]

, Trf = 1, (4.2)

where f is a generalized probability density. Equation (4.2) equals the grand canonicalpotential in equilibrium f = f0, with

Ω [f0] = −β−1lnΞ ≡ Ω, Ω [f ] > Ω [f0] , if f 6= f0. (4.3)

A proof of Eq. (4.3) is sketched in Appx. B.1.1.We consider a Hamiltonian of the form

HN

(

rN ,pN ,JN ,ωN)

=Ekin

(

pN ,JN)

+ U(

rN ,ωN)

+ V(

rN ,ωN)

(4.4)

=Ekin

(

pN ,JN)

+ U(

rN ,ωN)

+N∑

i=1

Vext (ri,ωi) , (4.5)

19

4. Density functional theory

where we decoupled the kinetic and potential energy analogous to Secs. 2.1 and 2.4. Ekin isthe kinetic contribution, U gives the potential energy resulting from the particle interactions,V represents the contribution of the external potentials and Vext is the external potentialacting on each individual particle i. Under the assumption of a system embedded in avolume we can define an equilibrium density as follows

ρ0 (r,ω) = 〈ρ (r,ω)〉 =

⟨

N∑

i=1

δ (r − ri) δ (ω − ωi)

⟩

= Trgc

[

f0

N∑

i=1

δ (r − ri) δ (ω − ωi)

]

, (4.6)

where ρ (r,ω) is the density operator (depending also on rN and ωN ) and δ is the Diracδ-distribution. Going back to the external potential, which we can rewrite as follows

V(

rN ,ωN)

=

∫∫

d3rd3ωρ (r,ω) Vext (r,ω) , (4.7)

Using Eq. (4.3) we find that

δΩ

δVext (r,ω)=

−1

βΞ

δΞ

δVext (r,ω)= ρ0 (r,ω) . (4.8)

In Eq. (4.8) we use the functional derivation (indicated by δ) as explained in Appx. B.3.Now it follows, that ρ0 (r,ω) is a functional of Vext, because f0 is a functional of Vext

itself. Further, f0 is a functional of ρ0 (r,ω), the proof is given in Appx. B.1.2. With theknowledge so far and using Eq. (4.8) we find that

ΩV [ρ] =

∫∫

d3rd3ωρ (r,ω)Vext (r,ω) + F [ρ] − µ

∫∫

d3rd3ωρ (r,ω) , (4.9)

with the intrinsic free energy

F [ρ] = Tr[

f0

(

Ekin + U + β−1lnf0

)]

, (4.10)

which is a unique functional of ρ (r,ω) for a given U .At this point we find that

ΩV [ρ] > ΩV [ρ0] = Ω,δΩV [ρ]

δρ (r,ω)

∣

∣

∣

∣

ρ0(r,ω)

= 0, (4.11)

here δ denotes a functional variation (see Appx. B.1.3 for the proof). Equation (4.11) isthe final expression for our density functional theory derivation, where the functional isdetermined by the one-particle density.

More general, the m-particle density distributions are derived as follows

ρ(m) (rm,ωm) = Ξ−1∞∑

N≥m

eβµN

ΛN (N −m)!

∫

d3rm+1 . . . d3rN

∫

d3ωm+1 . . . d3ωN ·

· e−β(V (rN ,ωN)+U(rN ,ωN)), (4.12)

20

4.2. Ideal gas

where Λ results from the integration over momenta and angular momenta, for more de-tails see Sec. 4.2. As mentioned earlier, the integration over the Euler angles reduces to∏N

j=m+1 dφjdθjsinθj in the case of linear molecules on which we will focus here. Within the

hierarchy of the density distributions we find ρ(1) (r,ω) ≡ ρ0 (r,ω). Further, we will use aproduct ansatz [28] for the density distributions, which is

ρ (r,ω) = ρ (r)α (r,ω) ,

∫

d2ωα (r,ω) = 1. (4.13)

Here α (r,ω) is the normalized distribution of molecular orientations. With this productrepresentation we find that in equilibrium the variations with respect to both ρ (r) andα (r,ω) must hold, that is

δΩ [ρ, α]

δρ (r)

∣

∣

∣

ρ0(r)= 0 (4.14)

and

δΩ [ρ, α]

δα (r,ω)

∣

∣

∣

α0(r,ω)= λ (r) , (4.15)

where λ (r) is a Lagrange multiplier, which ensures Eq. (4.13).

4.2. Ideal gas

We start with the free energy Eq. (2.18) and assume linear molecules, which possess n = 5degrees of freedom, i.e. 3 translational and 2 rotational. We can perform the momentumintegration d3Np = ΠN

i=1d3pi and angular momentum integration d3NJ = ΠN

i=1d3Ji as well

as the integration over the angles d2Nω = ΠNi=1dφidθsinθi, because the Hamiltonian is

separable [see Eq. (2.33)]. It follows then that

ZN (T, V ) =(4π)N

N !√

βh2/(2πm)3N√

βh2/(2πI)2N

∫

d3Nre−βVN . (4.16)

The remaining integration d3Nr = ΠNi=1d

3ri yields

Z idN =

1

N !

(

4πV

Λ3

)N

, (4.17)

where we use the fact of missing interactions, i.e. V idN = 0, in the ideal gas (id) and we

define Λtrans =√

h2β/ (2mπ), Λrot =√

h2β/ (2Iπ) and Λ = Λtrans (Λrot)23 . After inserting

the partition function in the free energy Eq. (2.19) and applying the Stirling approximation,

ln (N !) ≈ N lnN −N (4.18)

for large N , we arrive at

F id = −Nβ−1(

ln(

4πρΛ3)

− 1)

. (4.19)

Allowing radial and angular dependency and taking F id/V into account, we generalize byintegrating over the volume and the angles

F id [ρ] = β−1∫

d3r

∫

d2ωρ (r,ω)(

ln(

4πρ (r,ω) Λ3)

− 1)

. (4.20)

21


4.3. Fundamental measure theory

The aim of this section is to derive a density functional for the excess free energy Fex [ρ] foran inhomogeneous classical fluid composed of impenetrable spheres (the excess free energydue to anisotropic interactions is later discussed in Sec. 4.4). Contrary to the ideal gasin Sec. 4.2 it is much more cumbersome to deal with the excluded volume effects of hardspheres. In fact, a huge variety of approximations has been made, which are mostly basedon the incorporation of the analytic solution of the Percus-Yevick result for the direct paircorrelation function [56],

c(1) (r) = −β δFex [ρ]

δρ (r). (4.21)

There also is a variety of weighted density approaches [56]. Here, we will derive the fun-damental measure theory (FMT) based on Rosenfeld [57] and follow the description ofRoth [59]. The fundamental measure theory for mixtures seems to give the most promisingdescription of hard spheres and is even generalized to other than spherical shapes, e.g. inRef. [67]. The fundamental measure theory was often modified in the attempt to dealwith problematic issues arising from the dimensional crossover, that is the reduction ofthe density functional to lower dimensions by factorizing the density. We will discuss thedimensional crossover at the end of Sec. 4.3.3.

The derivation of the FMT is given in the following subsections. We start in Sec. 4.3.1with the so-called deconvolution of the Mayer-function, which will lead us to the originalRosenfeld functional [57] in Sec. 4.3.3 after a short consideration of the underlying scaledparticle theory in Sec. 4.3.2. The more present description of the White-Bear functional[58] is given in Sec. 4.3.4.

4.3.1. Starting point of a fundamental measure theory

The following derivation essentially follows Ref. [59] for hard-sphere mixtures, but we willlimit ourselves to the one-component system in this work. First we consider the grandcanonical potential,

−βΩ = lnΞ = ln

( ∞∑

N=0

1

N !Λ3N

∫∫

d3Nrd2Nωe−β 1

2

∑N

i=1

∑N

j=1,j 6=iφ(rij ,ωi,ωj)

· e−β∑N

i=1φext(ri,ωi)eβµN

)

, (4.22)

where we have already performed the integrations over momenta and angular momenta[cf. Eq. (4.2)]. In the next step we consider the Taylor expansion of

ln(1 +X) =∞∑

n=1

(−1)n−1Xn

n, (4.23)

which we use to rewrite Eq. (4.22) in orders of

z(ri,ωi) = eβµe−βφext(ri,ωi) = eβu(ri,ωi). (4.24)

22


This yields

−βΩ =∞∑

l=1

Λ−3lbl, (4.25)

where the first bl are

b1 =

∫∫

d3rd2ωz(r,ω), (4.26)

b2 =1

2

∫∫∫∫

d3rid3rjd

2ωid2ωjz(ri,ωi)z(rj ,ωj)f (rij,ωi,ωj) , (4.27)

b3 =1

6

∫∫∫

d3rid3rjd

3rk

∫∫∫

d2ωid2ωjd

2ωkz(ri,ωi)z(rj ,ωj)z(rk,ωk)

· f (rij ,ωi,ωj) f (rik,ωi,ωk) f (rjk,ωj ,ωk)

+1

2

∫∫∫

d3rid3rjd

3rk

∫∫∫

d2ωid2ωjd

2ωkz(ri,ωi)z(rj ,ωj)z(rk,ωk)

· f (rij ,ωi,ωj) f (rjk,ωj ,ωk) . (4.28)

At this point we introduce the Mayer-function,

f (rij,ωi,ωj) = e−βφ(rij ,ωi,ωj) − 1. (4.29)

In fact, Eq. (4.25) is still not very useful, since we are interested in an expression in ordersof the density ρ (ri,ωi). Therefore, one has to derive

ρ (ri,ωi) = − δΩ[u]

δu (ri,ωi)= −

∫∫

d3rd2ωδΩ[z]

δz (r,ω)

δz (r,ω)

δu (ri,ωi)= −z (ri,ωi)

δβΩ[z]

δz (ri,ωi).

(4.30)

After inserting Eq. (4.25) in Eq. (4.30) one has to solve the resulting equation with respectto z (ri,ωi). Motivated by the ideal gas solution zid = ρidΛ3, which results from

0 =d(βΩ/V )

dρ

∣

∣

∣

∣

ρid

= ln(

ρidΛ3)

− βµ, (4.31)

we assume z = ρΛ3 + correction. For the inhomogeneous anisotropic case we make theassumption

z(ri,ωi) =∞∑

l=1

(ρ(ri,ωi))l al, (4.32)

which is inserted in Eq. (4.30). Then comparing the coefficients leads us to the solution ofEq. (4.32). At this point we can insert Eq. (4.32) in Eq. (4.25) and rearrange the result inorders of ρ(ri,ωi). This leads us to the expansion of the excess (with respect to the idealgas) free energy at low densities,

βFex [ρ] = − 1

2

∫∫∫∫

d3r1d3r2d

2ω1d2ω2ρ (r1,ω1) ρ (r2,ω2) f (r12)

− 1

6

∫∫∫

d3r1d3r2d

3r3

∫∫∫

d2ω1d2ω2d

2ω3ρ (r1,ω1) ·

· ρ (r2,ω2) ρ (r3,ω3) f (r12) f (r13) f (r23) + O(

ρ4)

. (4.33)

23


In Eq. (4.33) we already use that we are interested in the hard-sphere pair potential, whichis independent of the orientation of the particles. The pair potential of hard spheres withdiameter σ = 2R is simply

φ (r) =

∞, |r| < σ,

0, |r| ≥ σ(4.34)

and therefore the Mayer-function reduces to

f (r) =

−1, |r| < σ,

0, |r| ≥ σ. (4.35)

As we see in the equation above, in our case the Mayer function is connected to the Heavisidefunction Θ (2R− r), for which we can find an expression by deconvolution, which is

Θ (2R− r) = −f (r) = 2 (ω3 ⊗ ω0 + ω2 ⊗ ω1 − ω2 ⊗ ω1) . (4.36)

The convolution of these scalar ωα and vector weight functions ωα is defined as

ωα ⊗ ωβ =

∫

d3r′ωα(

r′ − r)

ωβ

(

r′) , (4.37)

ωα ⊗ ωβ =

∫

d3r′ωα(

r′ − r)

ωβ

(

r′) . (4.38)

The weight functions which appear in Eq. (4.36) are

ω3 (r) = Θ (R− |r|) , (4.39)

ω2 (r) = δ (R− |r|) , (4.40)

ω1 (r) =ω2 (r)

4πR, (4.41)

ω0 (r) =ω2 (r)

4πR2, (4.42)

ω2 (r) = δ (R− |r|) r

|r| , (4.43)

ω1 (r) =ω2 (r)

4πR. (4.44)

As a side note: To get more used to this deconvolution and the weight functions we considerthe integrated Eq. (4.36),

−∫

d3rf (r) = 2 (V + SR) , (4.45)

which is just the excluded volume of two joined convex bodies. The integration over theweight functions ω3, ω2, ω1 and ω0 gives the volume, surface area, radius and the constant1, respectively, whereas the integral over the vector weight functions vanishes. These funda-mental geometric values give fundamental measure theory its name. Based on these weightfunctions we can define the so-called weighted densities

(

nα (r)

nα (r)

)

=

∫

d3r′ρ(

r − r′)(

ωα (r′)

ωα (r′)

)

, (4.46)

24


with which we can reformulate the first contribution to the excess free energy of Eq. (4.33)as follows

limρ→0

βFex [ρ, α] = − 1

2

∫∫∫∫

d3r1d3r2d

2ω1d2ω2ρ (r1,ω1) ρ (r2,ω2) f (r12)

=

∫

d3r (n0 (r)n3 (r) + n1 (r)n2 (r) − n1 (r) n2 (r)) . (4.47)

This ansatz so far does not reproduce higher-order terms in the excess free energy. There-fore, Rosenfeld suggested a more general expression, which gives us a functional for aninhomogeneous hard-sphere fluid of the form

Ω [ρ, α] =F id [ρ, α] + β−1∫

d3rφFMT (nα (r))

+

∫∫

d3rd2ωρ (r)α (r,ω)φext (r,ω) − µ

∫∫

d3rd2ωρ (r)α (r,ω) , (4.48)

using Eq. (4.13). To derive an expression for φFMT we perform a dimensional analysis, i.e.we require all terms in φFMT to be of dimension [length]−3. Thus Rosenfeld [57] formulatedthe expression

φFMT (nα (r)) =f1 (n3)n0 + f2 (n3)n1n2 + f3 (n3) n1n2 + f4 (n3)n32

+ f5 (n3)n2n22. (4.49)

It is worth mentioning, that another deconvolution of the Mayer function [68] in absenceof vector weight functions exists. A comparison of these two different deconvolutions givesus equations for f3 and f5 and yields

φFMT (nα (r)) = f1 (n3)n0 + f2 (n3) (n1n2 − n1n2) + f4 (n3)(

n32 − 3n2n2

2

)

. (4.50)

In the limit of vanishing density f1, f2 and f4 have to enforce the original deconvolution ofthe Mayer function. The next step will need the results from the so-called scaled particletheory, which is derived in the next section, Sec. 4.3.2.

4.3.2. Scaled particle theory

The original formulation of the FMT by Rosenfeld [57] incorporates a relation from thescaled particle theory to determine the functions fi (n3) in Eq. (4.50). Therefore, we de-scribe this theory in more detail.

The scaled particle theory, as derived in Ref. [69], gives us an approach for the work Wrequired to create a spherical cavity of radius R0 in a hard-sphere fluid. This work is thesame, as is required for inserting a solute sphere. We follow the steps of Ref. [69], wherethe derivation starts with the probability of spontaneous fluctuations forming such a cavity,which is

P0 (R0) = e−βW (R0). (4.51)

25


The probability that a sphere is in the cavity (overlap is prohibited) is

P1 (R0) = 1 − P0 (R0) =4π

3ρ (R0 +R)3 , (4.52)

where R is the diameter of a hard sphere, ρ is the number density of the fluid. Here R0

is used to vary the accessible volume for the fluid, i.e. 0 ≤ R0 + R ≤ R. For R0 ≤ 0 wecombine the equations above and find the following expression for the work:

W (R0) = −β−1ln

(

1 − 4πρ

3(R0 +R)3

)

. (4.53)

For R0 ≫ 0 we know, that the work of creating a huge cavity is given by the reduction ofthe accessible volume by ∆V0 under the pressure P , i.e.

W (R0) = P∆V0 =4π

3PR3

0. (4.54)

Assuming that we can write a cubic polynomial for the work if R0 ≥ 0, we incorporate thelimit case of Eq. (4.54) and get

W (R0) = w0 + w1R0 +w2R

20

2+

4πPR30

3. (4.55)

The coefficients ωi (i = 0, 1, 2) can be derived using the condition that Eq. (4.53) andEq. (4.55) are continuous in both the value and the first derivative by R0 = 0,

βw0 = − ln (1 − η) , (4.56)

βw1 =4πρR2

1 − η, (4.57)

βw2 =8πρR

1 − η+

(

4πρR2)2

(1 − η)2 , (4.58)

where the packing fraction η is 4πρR3/3. Now we can write down the excess chemicalpotential, µex, which equals the work for inserting a hard sphere of radius R = R0 into thefluid, the equation reads:

βµex = βW (R0) = −ln (1 − η) +6η

1 − η+

9η2

2 (1 − η)2 +βPη

ρ. (4.59)

The approximations appear on the right hand side of Eq. (4.59). Considering the limit ofinserting an infinitely large sphere we find

limR→∞

βµex

V=βP. (4.60)

26


4.3.3. The original Rosenfeld functional

Rosenfeld [57] used the expression for the excess chemical potential, Eq. (4.60), of the scaledparticle theory given in Sec. 4.3.2 to determine fi (n3) in Eq. (4.50). We will retrace thesteps of Ref. [59]. In fact, µex can be derived for a constant number density ρ in the absenceof an external field in the isotropic system [α (r,ω) = 1/ (4π)] by determining the derivativeof Eq. (4.48):

∂Ωbulk/V

∂ρ

∣

∣

∣

∣

ρ=ρ0

= − βµex +∂φFMT

∂ρ= 0. (4.61)

Reformulating this equation yields

βµex =∂φFMT

∂ρ=∑

α

∂φFMT

∂nα

∂nα

∂ρ, (4.62)

where we find that in the limit R → ∞ the derivative with respect to n3 ∝ R3 dominatesthe behavior, see Eqs. (B.38) and (B.39). Bringing all together and taking Eq. (4.60) intoaccount it follows that

limR→∞

βµex

V=∂φFMT

∂n3= βP. (4.63)

By substituting the thermodynamic bulk relation Eq. (2.24) for Eq. (4.48) in absence of anexternal potential, we find

−P =Ωbulk

V= φFMT +

F id

V− ρµ, (4.64)

which yields

βP = − φFMT +∑

α

∂φFMT

∂nαnα + n0. (4.65)

Connecting this result with Eq. (4.63) gives us a scaled particle differential equation for ourhard spheres, which is

∂φFMT

∂n3= − φFMT +

∑

α

∂φFMT

∂nαnα + n0. (4.66)

We solve this differential equation by collecting all contributions to each nα, integratingthe arising equations for every nα and choosing the integration constants such that the lowdensity limit, Eq. (4.47), holds. This gives us the original Rosenfeld functional,

φFMT (nα (r)) = − n0 (r) ln (1 − n3 (r)) +n1 (r)n2 (r) − n1 (r) n2 (r)

1 − n3 (r)

+n2

3 (r) − 3n2 (r) n22 (r)

24π (1 − n3 (r))2 . (4.67)

27


With this functional and Eq. (4.65) we can derive the Percus-Yevick equation of state

βPPY =n0

1 − n3+

n1n2

(1 − n3)2 +n3

2

12π (1 − n3)3 . (4.68)

Alternatively, Eq. (4.68) can be obtained by a functional variation of the grand canonicalfunctional [cf. Eq. (4.48)] with respect to ρ. The resulting expression is reformulated todetermine the chemical potential. The next step is to insert the chemical potential in thefunctional Eq. (4.48) and to use Eq. (2.24).

This original Rosenfeld functional describes the fluid phase very accurately, but failsat the dimensional crossover. In fact, strong confinement can be understood in the zero-dimensional limit, where a particle is trapped in a small cavity, therefore failure in thedescription of the solid phase is not surprising [59]. Addressing this, we give a brief dis-cussion of the dimensional crossover. As mentioned earlier, the density profile is factorizedsuch that it depends on the spatial coordinates of the lower dimension multiplied with adelta distribution depending on the higher dimension. Therefore, the assumption is that ifa functional is exact, it will be reduced to the exact functional in lower dimensions. Butin the one-dimensional limit the third term of the original Rosenfeld functional [57] [seeEq. (4.67)] diverges. This problem is then solved by considering the zero-dimensional limit,where only a small cavity is taken into account. Here the result is known from thermody-namics. The functional in higher dimensions can be modified to reproduce this result andthe problem arising in one-dimension is solved [70]. This modification of the functionalcauses different bulk results in three dimensions. Therefore, antisymmetrization can beused [70] [see Eq. (4.69)], which results in only small deviations in the zero-dimensionallimit. This antisymmetrization of the original functional leads to a functional which candescribe the solid phase, it reads as follows

φasy (nα (r)) = − n0 (r) ln (1 − n3 (r)) +n1 (r)n2 (r) − n1 (r) n2 (r)

1 − n3 (r)

+

n23 (r)

(

1 −∣

∣

∣

∣

n2 (r)

n2 (r)

∣

∣

∣

∣

)3

24π (1 − n3 (r))2 . (4.69)

Another way to solve these problems is the incorporation of tensorial weight functions asdone by Tarazona [71].

4.3.4. White-Bear version mark II

We choose the White-Bear version mark II functional [58] from the huge variety of availablehard-sphere descriptions in the framework of FMT. Hansen-Goos and Roth [58] start thederivation of the functional based on the Carnahan-Starling-Boublík equation of state

βPCS =n0

1 − n3+

n1n2

(1 − n3)2 +n3

2

12π (1 − n3)3 − n3n32

36π (1 − n3)3 . (4.70)

In contrast, Rosenfeld [57] focuses on the Percus Yevick equation of state [see Eq. (4.68)].Equation (4.65) defines a differential equation, where arising integration constants are fixed

28

4.4. DFT approximations to handle pair interactions beyond a reference potential

using the low density limit and the third viral coefficient which leads to the White-Bearfunctional (see Refs. [59] and [58]). This precursor version of the later mark II functional isquite successful and capable of describing crystals, but the exact thermodynamic relation[Eq. (4.60)] does not hold for the White-bear version [72]. Therefore, as described inRef. [59], an iterative loop is used to minimize the deviation from this relation. We omitthe details here and focus on the result. It follows that the numerical equation of state,

βPCSIII =n0

1 − n3+n1n2

(

1 + 13n

23

)

(1 − n3)2 +n3

2

(

1 − 23n3 + 1

3n23

)

12π (1 − n3)3 , (4.71)

minimizes this inconsistency, if it is used instead of the differential equation Eq. (4.65).Further, Eq. (4.60) then results in the Carnahan-Starling-Boublík equation of state. Solvingthis differential equation, the White-Bear functional mark II arises, which is

φWBII (nα (r)) = − n0 (r) ln (1 − n3 (r))

+n1 (r)n2 (r) − n1 (r) n2 (r)

1 − n3 (r)

(

1 +ψ2 (n3 (r))

3

)

+n2

3 (r) − 3n2 (r) n22 (r)

24π (1 − n3 (r))2

(

1 − ψ3 (n3 (r))

3

)

, (4.72)

ψ2 (n3 (r)) =

(

2n3 (r) − n32 (r) + 2 (1 − n3 (r)) ln (1 − n3 (r))

)

n3 (r), (4.73)

ψ3 (n3 (r)) =1

n32 (r)

(

2n3 (r) − 3n32 (r) + 2n3

3 (r)

+2 (1 − n3 (r))2 ln (1 − n3 (r)))

. (4.74)

It should be stressed here that the original derivation of all these functionals focuses onmixtures of two to three components. Therefore, a generalized Carnahan-Starling-Boublíkequation of state is used in Ref. [58].

4.4. DFT approximations to handle pair interactions beyond a

reference potential

Our attempts so far limit us to hard spheres. In a more complex system we have toincorporate free energy contributions due to additional particle interactions. This is adifficult task and therefore more approximations are necessary. From now on we assume apairwise additive particle interaction potential

U(

rN ,ωN)

=1

2

∑

i,j 6=i

φ (ri, rj ,ωi,ωj) . (4.75)

With such a pair interaction potential we find an expression for the two body densitydistribution which is different from Eq. (4.12) with m = 2:

2δΩ

δφ (ri, rj ,ωi,ωj)=

−2

βΞ

δΞ

δφ (ri, rj ,ωi,ωj)= ρ

(2)0 (ri, rj ,ωi,ωj) . (4.76)

29


We follow the Refs. [55] and [69] to define the so-called λ-expansion for the excess freeenergy functional

F [ρ0] = Fref [ρ0] +1

2

∫ 1

0dλ

∫∫∫∫

d3rid3rjd

2ωid2ωjρ

(2) (φλ, ri, rj ,ωi,ωj)

· ∂φλ (ri, rj ,ωi,ωj)

∂λ, (4.77)

which is obtained from functional integration of Eq. (4.76). The integration path starts ata reference potential and finishes at the complete pair interaction potential. We introduceλ = 0, . . . , 1 to describe this integration path and the pair potential for each value λ, whichis

φλ (ri, rj ,ωi,ωj) = φref (ri, rj) + λ(

φ (ri, rj ,ωi,ωj) − φref (ri, rj))

. (4.78)

In the equation above we assume that the reference potential does not depend on theorientations. This condition is satisfied for the hard-sphere potential, which is our choicefor the reference system.

4.4.1. Mean-field approximation

At this point we introduce an approximation for the treatment of the complex pair inter-actions shown above, Sec. 4.4. First, we use the definition of the pair correlation functiong (φλ, ri, rj ,ωi,ωj) to obtain an expression for the two particle density distribution,

ρ(2) (φλ, ri, rj ,ωi,ωj) = g (φλ, ri, rj ,ωi,ωj) ρ (ri,ωi) ρ (rj ,ωj) . (4.79)

In the equation above we can see that all correlational effects are included in this paircorrelation function. Therefore, the simplest approximation is to omit any correlationsexceeding those of the reference system. Often this procedure is described as setting g = 1,or in a more accurate notation

g (φλ, ri, rj ,ωi,ωj) ≈ g(

φref, ri, rj

)

. (4.80)

The reason for such an approach is the hope that some integrations decouple and the lackof an explicit analytic expression for the complete pair correlation function. Assuming thehard-sphere pair interaction as the reference system we find that the reference correlationfunction is just the Heaviside function Θ (|ri − rj | − σ), which is 1 for |ri − rj | ≥ σ and 0otherwise. Then only coupling between positions is gained from the pair potential, whichresults in the possibility to solve one integration over particle orientations in our model(see Sec. 3). Here we treat the anisotropic pair interaction in the mean-field approach anddescribe the reference system through the White-Bear mark II representation of a hard-sphere system [see Eq. (4.72)]. In mean-field approximation the resulting grand canonical

30

4.5. Confinement and walls

functional can be written as

Ω [ρ, α] =β−1∫

d3rρ (r)(

ln(

ρ (r) Λ3)

− 1)

+ β−1∫∫

d3rd2ωρ (r)α (r,ω) ln (4πα (r,ω))

+

∫∫

d3rd2ωρ (r)α (r,ω)φsurf (r,ω) − µ

∫

d3rρ (r)

+ β−1∫

d3rφWBII (nα (r))

+1

2

∫∫∫∫

d3rd3r′d2ωd2ω′ρ (r) ρ(

r′)α (r,ω)α(

r′,ω′)

·(

φ(

r1, r2,ω,ω′)

− φHS (r1, r2))

. (4.81)

4.4.2. Modified mean-field approximation

The modified mean-field theory (MMFT) and its application can be found in Refs. [73], [74]and [75]. This approach is found by considering the Ornstein-Zernike relation [69], fromwhich it is possible to derive that the pair correlation function in the low density limit is aBoltzmann factor. This yields

ρ(2) (φλ; ri, rj ,ωi,ωj) = g (φλ, ri, rj ,ωi,ωj) ρ (ri,ωi) ρ (rj,ωj) (4.82)

≈ ρ (ri,ωi) ρ (rj ,ωj) e−βφλ(ri,rj ,ωi,ωj). (4.83)

A more sophisticated approach incorporates an additional product of an exponential func-tion with further correlation functions. The reader may take a closer look in Ref. [69] formore details. As we will explain in Sec. 4.6, this approach gives rise to numerical computa-tion time. Therefore, we only apply it to the homogeneous isotropic case (see Sec. 4.6.1).

A comparison with Sec. 4.4.1 shows that the mean-field approach can be understood asthe high temperature limit of the modified mean-field approximation.

4.5. Confinement and walls

The incorporation of walls or the description of confined geometry enforces inhomogeneousdensity profiles. If, however, one wall is located at z = 0 the density profile in the limit ofinfinite distance from the wall approaches the constant bulk density of a fluid, which means

ρ (z → ∞) = ρbulk. (4.84)

Taking a closer look at the behavior close to the wall, we consider the contact theoremfor planar walls [56]. This theorem connects the pressure P of the fluid to the force perunit area caused by the wall on the fluid.

P = −∫∫

dzd2ωρ (r · n,ω) F (r,ω) · n

= −∫

dzρ (z)

∫

d2ωα (z,ω)∂

∂zφsurf (z, u (ω)) , (4.85)

31


where n is the normal unit vector on the surface, here in z-direction. Further, F (r,ω)is the force exerted by the wall on the fluid, and φsurf (z,ω) is the surface potential. Byconsidering a hard wall at z = 0 the contact theorem is reduced to [69]

P = β−1∫ ∞

0dz

d

dz

(

e−βφsurf)

y (z) = β−1∫ ∞

0dz

d

dz

(

Θ

(

z − σ

2

))

y (z) . (4.86)

Here, we assume a density function, such as

ρ (z) = e−βφsurfy (z) , (4.87)

where y (z) is a continuous function in z [69]. The contact theorem for a hard planar wallthen is

ρ

(

z =σ+

2

)

= βP bulk, (4.88)

where the + indicates z = ǫ + σ/2 with ǫ infinitely small. This theorem expresses thebalance of forces, so the force per unit area exerted by the wall equals the bulk pressure.

4.6. Density functional theory applied to our model

Our starting point is Eq. (4.81) in the mean-field approach, where we use the specificproperties of Eq. (3.5) as shown in the Appx. B.4.1, yielding

Ω [ρ, α] =β−1∫

d3rρ (r)(

ln(

ρ (r) Λ3)

− 1)

+ β−1∫∫

d3rd2ωρ (r)α (r,ω) ln (4πα (r,ω))

+

∫∫

d3rd2ωρ (r)α (r,ω) φsurf (r, u (ω))

− µ

∫

d3rρ (r) + β−1∫

d3rφWBII (nα (r))

+

∫∫∫

d3rd3r′d2ωρ (r) ρ(

r′)

· α (r,ω)φ1(∣

∣r − r′∣

∣

) r′ − r

|r − r′| · u (ω) . (4.89)

Of course the whole formalism can be written in the modified mean-field approximation,but it would become necessary to expand the α (r,ω) in Legendre polynomials and rewritethe pair potential in spherical invariants again. Then the number of expansion coefficientsneeded to improve the approach compared to a simple mean-field approximation woulddrastically increase the computation time. One also needs to consider that the expansioncoefficients of α (r,ω) have to be extracted at the end, which will give rise to two additionalintegrations over the angles. We limit ourselves to the study of the homogeneous isotropicsystem in the modified mean-field approach in Sec. 4.6.1, focusing on a possible liquid-vaporphase transition.

32

4.6. Density functional theory applied to our model

Going back to Eq. (4.89): Treating the anisotropic contribution to the functional inmean-field approximation, we are interested in the equilibrium solution, which minimizesthe functional. Therefore we calculate the functional derivative with respect to α (r,ω),which leads to

α0 (r,ω) =exp (a (r) · u (ω))

∫

d2ω′exp (a (r) · u (ω′)), (4.90)

where

a (r1) = β

∫

d3r2ρ (r2)φ1 (r12) r12 − βφsurf (r1 · ez) . (4.91)

The complete derivation can be found in the Appx. B.4.2. The next step is to reinsert the so-lution, Eq. (4.90), into the functional, Eq. (4.89), and to determine the functional variationwith respect to ρ (r). In our study we are interested in the competition between fluid-fluidand fluid-surface interactions, therefore we consider a planar surface located at z = 0. Werestrict ourselves to ρ (r) = ρ (z) due to considerations of numerical effort. Clearly, thisstrategy implies the possibility of missing other, maybe energetically preferable, solutions.Indeed, by considering both z-dependent and r-dependent density distributions, Tarazonaand coworkers [28] showed that the present model yields not only planar structures (“mem-branes”), but also spherical structures (“vesicle” and “micelles”) under bulk conditions.However, in the presence of surfaces a full three dimensional calculation would result inconsiderable computational cost. Therefore, we focus on density distributions with thesymmetry suggested by the planar wall and assume translational symmetry in the otherspatial dimensions. The functional derivation with respect to ρ (z) can be found in theAppx. B.4.3. Here we will consider the final result, which is

ρ (z) = eνρ0e−βφsurf,hard(z)−β 1

A

δFex[ρ]δρ(z) . (4.92)

In this equation Fex[ρ] is the excess free energy functional. It includes both the HS and theanisotropic interaction contributions to the full grand canonical density functional equation,as well as the contribution of the orientation dependent part of the external (surface) po-tential. Furthermore, A denotes the surface area, φsurf,hard (z) is the potential contributionfor z < σ/2, and νρ0 is defined in Eq. (B.43). This final expression has to be evaluatednumerically, our algorithm and some advice are given in the Appx. B.2.

4.6.1. The homogeneous isotropic system in modified mean-fieldapproximation

The aim of this subsection is to derive the pressure P = − Ω

V= PHS + P an for the homo-

geneous isotropic system, where P an is the pressure contribution caused by the anisotropicpair interaction. Therefore, we derive the equilibrium solution for the functional Ω [ρ]through functional variation. We gain an expression for the chemical potential, which weinsert in the original functional, which is then in equilibrium. Contrary to our inhomo-geneous density functional theory study we use the original Rosenfeld functional for the

33


hard-sphere contribution, which gives us the Percus-Yevick pressure Eq. (4.68). The mod-ified mean-field approach, Eq. (4.83), using Eq. (4.77) for the anisotropic contribution ofthe pair interaction, Eq. (3.8), gives us the equation

Ωan [ρ, α] = − 1

2β

∫∫∫∫

d3r1d3r2d

2ω1d2ω2ρ (r1) ρ (r2)α (r1,ω1) ·

· α (r2,ω2) e−βφHS(r12)(

e−βφI(r12,ω1,ω2) − 1)

. (4.93)

With our assumptions ρ (ri) = ρ = const and α (ri,ωi) = (4π)−1 it follows that

Ωan [ρ] = 2πV ρ2β−1∫ ∞

σdrr2

(

sinh (−βφ1 (r))

βφ1 (r)+ 1

)

. (4.94)

After calculating the contribution of Ωan [ρ] to the chemical potential via derivation withrespect to ρ and inserting the chemical potential we get

βP = βPPY + 2πρ2∫ ∞

σdrr2

(

sinh (−βφ1 (r))

βφ1 (r)+ 1

)

. (4.95)

The integral above is negative, therefore a minimum and a maximum, or a saddle pointcan be found in the pressure as a function of the density. The occurrence of such a saddlepoint or maximum and minimum gives us a hint of a possible gas-liquid phase transition,which can be derived using the Maxwell construction (see Appx. A.3).

4.6.2. Parameters for our calculations

The free parameters for our system are the bulk packing fraction

η =π

6ρσ3 =

π

6ρ∗, (4.96)

the reduced temperature

T ∗ =kBTσ

C(4.97)

and the coupling strength to a surface potential C2. Obviously the T ∗ is inversely propor-tional to the anisotropic coupling constant C (see Eq. (3.8)) and in the limit of infinitelyhigh T ∗ the pair potential Eq. (3.5) will reduce to the repulsion of hard spheres. T ∗ & 10is a good approximation for this limit as we will see. The other parameters, λ∗ and λ∗

2,(mentioned earlier) are fixed to the value 3, if not stated otherwise.

34

5Chapter 5.

Molecular dynamics simulations

This chapter is structured as follows: We start with the derivation of our basic Verletalgorithm in Sec. 5.1, followed by the introduction of the Berendsen thermostat in Sec. 5.2.The handling of periodic boundaries is described in Sec. 5.3, while Sec. 5.4 refers to thetreatment of the interaction range of the potentials. The important point of achievingequilibrium and its characterization is addressed in Sec. 5.5. We briefly consider finite sizeeffects in Sec. 5.6 and close the chapter with some numerical details in Sec. 5.7.

5.1. The Verlet algorithm

We choose the velocity Verlet algorithm given in Ref. [76] to solve the equations of motion.The derivation starts with the Taylor expansion of the position ri (t+ ∆t) of particle i,which is

ri (t+ ∆t) = ri (t) +

(

ri (t) +∆t

2ri (t)

)

∆t+O(

∆t3)

, (5.1)

ri (t+ ∆t) = ri (t) + ri

(

t+∆t

2

)

∆t+O(

∆t3)

. (5.2)

We rename the velocity ri (t) as the vi (t), and we use the force Fi acting on particle i todetermine the acceleration ri (t) = Fi (t) /mi.

Rotational dynamics are incorporated by considering the Taylor expansion of the vectorof orientation ui (t+ ∆t),

ui (t+ ∆t) = ui (t) +

(

ui (t) +∆t

2ui (t)

)

∆t+O(

∆t3)

, (5.3)

ui (t+ ∆t) = ui (t) + ui

(

t+∆t

2

)

∆t+O(

∆t3)

. (5.4)

In general we would need to use a rotation matrix R to rotate the unit vector ui (t) =ui (t) / |ui (t)| around the rotational axis ui (t) × ui (t+ ∆t) yielding ui (t+ ∆t) = Rui (t).By allowing rotations through only small angles, we can approximate the rotation matrixand then express ui (t) as ωi (t) × ui (t) [77], where ωi (t) is the angular velocity. Thenit follows that ui (t) = ωi (t) × ui (t) [77], where we can derive the ωi (t) from the torqueT (t) = Iiωi (t) acting on a spherical particle i with the moment of inertia Ii.

By combining all we find the Verlet algorithm for translational and rotational degreesof freedom. The algorithm starts with the calculation of the new velocities and angular

35

5. Molecular dynamics simulations

velocities at a half-step, i.e.,

vi

(

t+∆t

2

)

= vi (t) +Fi (t)

mi

∆t

2, (5.5)

ωi

(

t+∆t

2

)

= ωi (t) +Ti (t)

Ii

∆t

2. (5.6)

In fact, the difference between vi (∆t/2) and vi (0) as well as the difference in the analogousangular velocities is often ignored [76], a practice which we choose to follow. The secondstep of the Verlet algorithm calculates the positions and the orientations, which are

ri (t+ ∆t) = ri (t) + vi

(

t+∆t

2

)

∆t, (5.7)

ui (t+ ∆t) = ui (t) + ωi

(

t+∆t

2

)

× ui (t) ∆t, (5.8)

ui (t+ ∆t) =ui (t+ ∆t)

|ui (t+ ∆t)| . (5.9)

Here we have to normalize the orientations at the end of this step [e.g., as done in thepopular program code LAMMPS (see Appx. D)]. This is needed because the approximationof the small angular rotation does not enforce the normalization by itself. The force andthe torque acting on a particle i follow as

Fi (t+ ∆t) =∑

j,j 6=i

Fij (t+ ∆t) = −∑

j,j 6=i

∇riφ (rij , ui, uj) , (5.10)

Ti (t+ ∆t) =∑

j,j 6=i

Tij (t+ ∆t) = −∑

j,j 6=i

ui (t+ ∆t) × ∇uiφ (rij, ui, uj) . (5.11)

The resulting expressions for the forces and torques of our model can be found in theAppx. C.2. After determining Eqs. 5.10 and 5.11 the next step is the final calculation ofthe velocities and angular velocities, i.e.,

vi (t+ ∆t) = vi

(

t+∆t

2

)

+Fi (t+ ∆t)

mi

∆t

2, (5.12)

ωi (t+ ∆t) = ωi

(

t+∆t

2

)

+Ti (t+ ∆t)

Ii

∆t

2. (5.13)

5.2. Thermostat

The native ensemble of molecular dynamics simulations is the microcanonical ensemble.Through the incorporation of a thermostat we can adapt the simulations to other ensemblesranging from the isokinetic to the canonical ensemble. A thermostat, which can easily beimplemented, is the Berendsen thermostat [78]. This thermostat couples the system to anexternal heat bath with constant temperature T0 and preserves the Maxwellian shape ofthe velocity profiles [78]. We follow the arguments of Ref. [78] and recapitulate the basicprinciples of this thermostat. The derivation of this thermostat is based on the Langevinequation,

mivi (t) = Fi (t) −miγvi (t) + Ri (t) , (5.14)

36

5.2. Thermostat

where vi is the time derivative of the velocity, γi is a friction constant and Ri represents aGaussian stochastic variable with

〈Ri (t)〉 = 0 (5.15)

and

〈Ri (t) Rj (t+ τ)〉 = 6miγikBT0δ (τ) δij . (5.16)

In the equation above δij is the Kronecker delta. It is assumed that γi = γ. The Langevinequation, Eq. (5.14), couples locally and globally to the external heat bath. Berendsen etal. focus on a global coupling to a heat bath, which is achieved by modifying Eq. (5.14).At first, Berendsen et al. [78] calculate the time-derivative of the kinetic energy as follows

dEkin

dt= lim

∆t→0

[

1

2∆t

N∑

i=1

miv2i (t+ ∆t) − 1

2∆t

3N∑

i=1

miv2i (t)

]

. (5.17)

Using the relation for the velocity change in one time step,

vi (t+ ∆t) = ∆vi + vi (t) =

∫ t+∆t

tdt′vi

(

t′)

+ vi (t) , (5.18)

and using the properties of the Gaussian stochastic variables, it follows that

dEkin

dt=

N∑

i=1

vi (t) · Fi (t) + 2γ

(

3N

2kBT0 − Ekin

)

(5.19)

=N∑

i=1

vi (t) · Fi (t) +Ekin

τT

(

T0

T− 1

)

(5.20)

with τ = (2γ)−1. The global coupling to a heat bath is achieved using the expressionfor the acceleration resulting from Eq. (5.20) of a particle i. It additionally utilizes thecorresponding modified Langevin equation [78],

mivi (t) = Fi (t) +mi

2τT

(

T0

T− 1

)

vi (t) , (5.21)

which can be expressed with respect to a time step:

vi (t+ ∆t) =Fi (t)

mi∆t+

(

1 +∆t

2τT

(

T0

T− 1

))

vi (t) . (5.22)

The equation above describes a rescaling of the velocities in every time step. We find thetemperature change per time step by inserting this result into the kinetic (here transla-tional) energy. By changing the rescaling factor Berendsen et al. [78] set this change to(T0 − T ) ∆t/τT . For the numerical implementation the following equation then holds

vi (t+ ∆t) = vi (t+ ∆t)

√

1 +∆t

τT

(

T0

T− 1

)

. (5.23)

37


The created ensemble is not rigorously defined [78]. Interestingly the authors write: “Wedo not consider this as a disadvantage because in cases where coupling to an externalbath is required, ..., ensemble fluctuations are generally insufficiently accurate to be usefulas a source for the deviation of macroscopic properties” [78]. We note that Morishitapublished a study [79] on fluctuation formulae for the Berendsen thermostat, which takeinto account that this thermostat does not reproduce an canonical ensemble. This influencesthe calculation of the specific heat capacity, see Sec. 6.5.

In our molecular dynamics simulations we also deal with rotations of our particles. Weconsider the torque

Ti (t) = Iiωi (t) = ui (t) × vi (t)mi, (5.24)

where ui is the orientation vector of a sphere with length of the radius and Ii representsthe moment of inertia for a spherical particle. At this point we see that the formalism ofa Berendsen thermostat, i.e. the rescaling of the velocities vi (for rotations the velocitytangential to ui), can also be applied to rotational dynamics by rescaling the angularvelocities ωi. For our simulation this yields

ωi (t+ ∆t) = ωi (t + ∆t)

√

1 +∆t

τT

(

T0

Trot− 1

)

. (5.25)

We can control the coupling strength of the system with the heat bath through the pa-rameter τT . In the limit of τT = ∆t the rescaling factor reduces to

√

T0/T an the createdensemble is isokinetic. Values of τ higher than the time step allow temperature fluctuations.We use τT = 10∆t. The temperature of the system is calculated using the kinetic energy.In equilibrium, the temperature calculated from the total kinetic energy equals each ofboth temperatures Ttrans and Trot, which are determined from the translational and therotational contribution to the kinetic energy. These two temperatures are given as

kBTtrans =

∑

i miv2i

(

t+∆t

2

)

3N − 4, (5.26)

kBTrot =

∑

i Iiω2i

(

t+∆t

2

)

2N − 1. (5.27)

The subtraction of degrees of freedom in the denominator of both formulae above takesinto account that, firstly, the total momentum is conserved and, secondly, the total kineticenergy is fixed.

There is a huge variety of thermostats available, e.g. the Anderson, the Nosé Hoover andthe constrained thermostat, the reader may take a closer look at Ref. [77]. In fact, theBerendsen thermostat is not a particularly sophisticated one, but easy to implement. Sincewe are mainly interested in equilibrium properties and structure analysis, the incorporationof other and maybe better thermostats is not considered to be crucially important.

38

5.3. Boundary conditions

0 1 2 3 4r/σ

10

20

30

0 1 2 3 4r/σ

0

2

4

6T*=1, rc

*=5

T*=0.1, rc

*=5

T*=1, rc

*=4

T*=0.1, rc

*=4

0 1 2 3 4r/σ

0

200

400

600

g 000(r

)

b)a) c)

Figure 5.1.: Radial pair correlation function at a) ρ∗ = 0.01, b) ρ∗ = 0.1, c) ρ∗ = 0.5 at tworeduced temperatures T ∗, a detailed discussion of the results follows in Sec. 7.2.

5.3. Boundary conditions

In the absence of any confining wall in a given spatial direction we apply periodic boundaryconditions, which connect exact copies of the simulation box in this spatial direction. Thenumerical implementation lets every particle leaving the box on one side enter the same boxat the opposite side. The drawback is that the total angular momentum is not conservedby this method.

5.4. Truncation and shifting of potentials

In general, one treats short and long range interactions in different ways. Consideringa truncation of the interaction at a cutoff radius rc, we omit a tail contribution of theinteraction potential, φ (r), to the average potential energy [80]

Etailpot = 2πN

∫ ∞

rc

drr2φ (r) ρ (r) . (5.28)

If we assume a constant density ρ (r) = ρ, Eq. (5.28), with φ (r) ∝ r−α, diverges for α ≤ 3.Hence, truncations should only be applied to short range interactions, where Etail

pot can beneglected. One definition of short range interactions is that they decay faster than r−6

[81]. Long range interactions on the other hand can be treated with the Ewald sum [77],which is probably one of the most sophisticated methods. In our model we deal with shortrange interactions, like the soft-sphere repulsion and the Yukawa potential (see Sec. 3.2),where we truncate the potential after a distance rc = 5σ. Considering our pair potentialfor various temperatures, we find the values at 5σ to be close to zero and at least 5 ordersof magnitude or more orders of magnitude smaller than the contribution at distance 1σ.This provides the numerical justification for this treatment. It should be noted, that theconsidered interaction length, the point of truncation, should not be larger than half thebox size, otherwise particles will interact with their periodic image.

39


0 5e+06 1e+07

t*

-20

-15

-10

-5

0

Eto

t* / N

Figure 5.2.: Total energy per particle as a function of the reduced simulation time t∗ atρ∗ = 0.01, T ∗ = 0.1.

This truncation causes a jump in the potential from a value close to zero to zero. Acommon way to omit the discontinuity in the potential and the resulting change in thetotal energy of the system is to shift the potential by its value at the truncation distance[77], when a particle approaches this boundary. Further, it may be necessary to shift theforces, too [77]. In fact, the influence of such deviations depends on the truncation pointand on the radial dependency of the pair potential. Here, we only shift the soft-spherepair potential and not the forces. This procedure is used in the popular program codeof LAMMPS (see Appx. D) for spherical particles with dipoles. In the following we showthe resulting total energies per particle, Etot, for some parameter combinations of theparameters T ∗ (reduced temperature), ρ∗ (reduced density) and r∗

c = rc/σ (see Appx. A.4).

T ∗ ρ∗ r∗c 〈Etot〉 T ∗ ρ∗ r∗

c 〈Etot〉 T ∗ ρ∗ r∗c 〈Etot〉

1 0.01 5 2.50 1 0.1 5 2.59 1 0.6 5 4.551 0.01 4 2.50 1 0.1 4 2.59 1 0.6 4 4.55

0.1 0.01 5 -20.10 0.1 0.1 5 -18.20 0.1 0.6 5 -15.430.1 0.01 4 -20.17 0.1 0.1 4 18.21 0.1 0.6 4 -15.45

We find only small deviations in the total energy, which are more pronounced for lowerreduced temperatures. The magnitude of the deviations is of the order of the standarddeviation of each mean value. These deviations are small enough to be ignored.

Further, we show the radial pair correlation function [see Eq. (6.5)] for the correspondingparameter pairs (T ∗, ρ∗, r∗

c ) in Fig. 5.1, which agree for both truncation lengths and all re-duced densities, ρ∗, considered. Summarizing, 4σ appears to be a sufficient large truncationpoint, we use the larger value anyway.

5.5. Characterizing equilibrium

The total energy in a simulation using a thermostat is not fixed, as it is in a microcanonicalsimulation, but in equilibrium the total energy fluctuates around a mean value. Therefore,

40

5.6. Finite size effects

0 5 10r/σ

0

0.5

1

1.5

2

g 000C

lust

er(r

)

0 2 4 6 8r/σ

0

1

2

3

4

5

6N = 512, T* = 0.15

N = 1000, T* = 0.15

N = 2197, T* = 0.15

N = 512, T* = 0.1

N = 1000, T* = 0.1

N = 2197, T* = 0.1

b)a)

Figure 5.3.: Finite size effects of the correlation function between two cluster centers at a)ρ∗ = 0.1 and b) ρ∗ = 0.5. A detailed discussion of the results follows in Sec. 7.2.

we can define a system as equilibrated, when the average value of the total energy overtime becomes constant. Additionally, time averages in equilibrium do not change (withinnumerical accuracy), so we also check the change in the radial distribution function withinthe simulation time. For most densities considered, equilibrium is found to set in within areasonable time. However, problems occur in very dilute systems, such as at the reduceddensity ρ∗ = ρσ3 = 0.01, and at low reduced temperatures (i.e. high anisotropic coupling).Such a case is shown in Fig. 5.2, where we plot the total energy as a function of thesimulation time for the low reduced temperature T ∗ = 0.1. Even after over 100 milliontime steps (∆t∗ = ∆t

√

ǫ/(σ2m) = 0.005) and after reaching temporary almost stable totalenergy plateaus the average of the total energy decreases further (though very slowly). Atthe same time, the first peak of the radial distribution function [see Eq. (6.5)] continues togrow, indicating slowly growing clusters. For such cases as just described the aggregationis not completely finished even after these very long simulation times. In all other cases,so at densities ρ∗ ≥ 0.1, these problems do not arise and we find equilibrium after severalmillion time steps. For safety we calculate at least 30 million time steps before starting anycalculation of time averages.

At low temperatures we additionally check the equipartition of the kinetic energy. Thekinetic energy contribution of each cluster is divided by the number of particles boundin the corresponding cluster. For ρ∗ ≥ 0.1 the overall clusters average of this quantityequals the kinetic energy contribution of a particle and we find the standard deviation tolie between 16% and 20% of the average value at T ∗ = 0.1.

5.6. Finite size effects

Due to restrictions on the computation time, the number of particles which can be incor-porated in the simulation is limited. Hence, finite size effects can occur in our simulations,which could be significant. A sophisticated finite size study would require a lot of computa-tion time, so due to the long equilibrating time (see Sec. 5.5) we concentrate on some testcalculations. Our simulation use 1000 particles, and we perform additional test simulations

41


0 5 10 15s

0

0.2

0.4

0.6

NC(s

)/ N

C

0 10 20 30s

0

0.2

0.4

0.6

0.8N = 512, T* = 0.15N = 1000, T* = 0.15N = 2197, T* = 0.15N = 512, T* = 0.1N = 1000, T* = 0.1N = 2197, T* = 0.1

a) b)

Figure 5.4.: Finite size effects of the fraction of clusters NC(s) as a function of cluster sizes and total number of clusters N tot

C for a) ρ∗ = 0.1, b) ρ∗ = 0.5, a detaileddiscussion of the results follows in Sec. 7.2.

with 512 and 2197 particles. We present some averages of the total energy per particle Etot

in the following table.

T ∗ ρ∗ N 〈Etot〉 T ∗ ρ∗ N 〈Etot〉0.17 0.01 512 -7.06 0.15 0.5 512 -4.710.17 0.01 1000 -3.13 0.15 0.5 1000 -4.710.17 0.01 2197 -1.96 0.15 0.5 2197 -4.720.1 0.01 512 -20.37 0.1 0.5 512 -16.740.1 0.01 1000 -20.10 0.1 0.5 1000 -16.960.1 0.01 2197 -18.13 0.1 0.5 2197 -16.960.19 0.1 512 -1.32 0.14 0.8 512 -1.430.19 0.1 1000 -1.31 0.14 0.8 1000 1.44)0.19 0.1 2197 -1.28 0.14 0.8 2197 -1.420.1 0.1 512 -18.27 0.1 0.8 512 -10.200.1 0.1 1000 -18.21 0.1 0.8 1000 -10.240.1 0.1 2197 -18.20 0.1 0.8 2197 -10.14

For most densities ρ∗ ≥ 0.1 the changes in the total energy per particle Etot are found tobe negligible. Surprisingly, for ρ∗ = 0.01 appreciable differences occur at low reduced tem-peratures, which may be related to the difficulties in achieving equilibrium at this density,as mentioned in Sec. 5.5. Considering the cluster center correlations [defined analogousto Eq. (6.5)] shown in Fig. 5.3 for two T ∗ and ρ∗ we find good agreement for all systemsizes N , with some derivations for ρ∗ = 0.5, T ∗ = 0.1 at large distances r. Further, weplot the cluster size distributions divided by the total number of clusters for two ρ∗ and T ∗

in Fig. 5.4. Again we find good agreement, but in the case of ρ∗ = 0.1 and T ∗ = 0.1 wesee deviations in the ratio of the most likely clusters of size 6 − 9. In the densest systemconsidered, ρ∗ = 0.8, we find deviations of the ratio between clusters of size s = 13 and 19(see Fig. 5.5a)) in the cluster size distribution. The clusters of size 13 dominate for N = 512.These differences also occur in the cluster center correlation shown in Fig. 5.5b), where werecognize a double peak structure around r = 3/σ for N = 1000 and only one peak shifted

42

5.7. Additional numerical details

0 5 10 15 20 25 30 35 40 45s

0

0.1

0.2

0.3

0.4

NC(s

)/ N

C

0 1 2 3 4 5 6 7r / σ

0

1.5

3

g 000C

lust

er(r

)

N = 512, T* = 0.1N = 1000, T* = 0.1N = 2197, T* = 0.1

a) b)

Figure 5.5.: Finite size effects of the a) fraction of clusters NC(s) as a function of clustersize s and total number of clusters N tot

C , b) correlation function of two clustercenters for ρ∗ = 0.8, a detailed discussion of the results follows in Sec. 7.2.

to the left for the other two system sizes. In fact, we cannot find a general increase ordecrease of the ratio between both cluster sizes s, while we are increasing the system sizeN . The reason is that the relative number of clusters of size s = 19 grows for N = 1000,but decreases again for N = 2197. We have to stress that for the cluster search (considering1 million time steps) we average with less accuracy for N = 512 and 2197. In Fig. 5.6 weshow the specific heat capacity as a function of T ∗ for two different ρ∗. Even though thefinite size effects in C∗

V (see Sec. 6.5) are not great, we see some deformation of the peakfor ρ∗ = 0.1 and a slight shift of the peak for ρ∗ = 0.5 to lower T ∗. As the simulations withN = 512 and 2197 run less long than those with N = 1000 some deviations in the errorsensitive C∗

V can maybe additionally result from the runtime difference.

A further important note is that we perform this finite size consideration for the “basic”model without extensions [see Eq. (3.12)] in the bulk system. It is possible that finite sizeeffects are stronger for the extended model Eq. (3.9) (q1 = 0), due to the formation of largeraggregates and even layers.

5.7. Additional numerical details

In our molecular dynamics simulations, the ratio of the thermal energy and the soft-sphererepulsion constant is set constant to kBT/ǫ = 1. This choice is made to focus on theinfluence of the anisotropic pair interaction potential, where the pure soft-sphere repulsionwould lead to a fluid phase for all densities considered (ρ∗ ≤ 0.8). Further, the ratio kBT/ǫ =1 allows only small interpenetration of spheres and maintains a qualitative comparabilityto hard spheres, as studied within our DFT formalism. To characterize the strength of theanisotropic interaction, we define a reduced temperature T ∗ = kBTσ/C.

Time averages over energies and related quantities are typically taken over a time intervalof 3 million time steps with 6000 equidistant and uncorrelated measurements. The timeaverages concerning cluster distributions and correlation functions are taken over 2 to 3million time steps. For time correlation functions we average over 100 time intervals with

43


0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3

T*

0

2

4

6

8

C* V

N = 512, ρ* = 0.1

N = 1000, ρ* = 0.1

N = 2197, ρ* = 0.1

N = 512, ρ* = 0.5

N = 1000, ρ* = 0.5

N = 2197, ρ* = 0.5

Figure 5.6.: Finite size effects of the specific heat capacity as a function of the reducedtemperature for two different reduced densities, a detailed discussion of theresults follows in Sec. 7.2.

different time origins.

44

6Chapter 6.

Correlation functions and furtherrelevant quantities

In this chapter we define and derive relevant quantities to analyze our simulation and DFTresults. We start with the structure correlation functions in Sec. 6.1, which are followed bythe definition of some additional order parameters and functions in Sec. 6.2. Further, weaddress time correlation functions in Sec. 6.3 and the mean-square displacement in Sec. 6.4.Finally, we introduce the specific heat capacity in Sec. 6.5 and the pressure in Sec. 6.6.

6.1. Structure correlation functions

The angular pair correlation function g (r,ω1,ω2) is our starting point for all structurecorrelation functions. We follow [60] and expand g (r,ω1,ω2) in rotational invariants withthe expansions coefficients gl1l2,l (r) for a homogeneous isotropic fluid. For linear moleculesthis expansion in the space fixed frame (coordinate system) reads

g (r,ω1,ω2) =∑

l1,l2,l

∑

m1,m2,m

gl1l2,l (r)C (l1, l2, l;m1,m2,m)Yl1,m1 (ω1)

· Yl2,m2 (ω2)Yl,m (ω)∗ , (6.1)

where C (l1, l2, l;m1,m2,m) are the Clebsch-Gordon coefficients, Yl,m (ω) is a sphericalharmonic and ∗ gives the complex conjugate. Then the projections can be determinedusing [60]

gl1l2,l (r) = 16π2∑

m

√

4π

2l + 1C (l1, l2, l;m,−m, 0)

·⟨

g(

r,ω′1,ω

′2

)

Yl1,m

(

ω′1

)∗Yl2,−m

(

ω′2

)∗⟩

ω′1,ω′

2

, (6.2)

where the prime denotes the intermolecular frame. The pair correlation function is definedas

g (r,ω1,ω2) =

(

4π

ρ

)2⟨

N∑

i,j 6=i

δ (ri) δ (rj − r) δ (ω1 − ωi) δ (ω2 − ωj)

⟩

(6.3)

= V

(

4π

N

)2⟨

N∑

i,j 6=i

δ (r − rij) δ (ω1 − ωi) δ (ω2 − ωj)

⟩

. (6.4)

The corresponding Clebsch-Gordon coefficients and the spherical harmonics are given inAppx. A.5. In the following, some constant factors (see Appx. A.5) will be omitted in the

45

6. Correlation functions and further relevant quantities

definition of our correlation functions. This choice leads to limr→∞ g000 (r) = 1 and thehigher-order correlation functions still vanish for r → ∞. The lowest order projections are

g000 (r) ≡⟨

V

4πr2N2

∑

i

∑

j 6=i

δ (r − rij)

⟩

, (6.5)

g110 (r) ≡⟨

V

4πr2N2

∑

i

∑

j 6=i

δ (r − rij) ui · uj

⟩

, (6.6)

g101 (r) ≡⟨

V

4πr2N2

∑

i

∑

j 6=i

δ (r − rij) ui · rij

⟩

. (6.7)

In the equations above δ (. . .) is the Dirac delta distribution and the term 4πr2 in thedenominator arises from δ (r − rij), where the restriction on the radial distance leads toδ (r − rij) /(4πr2). The radial distribution function g000 (r) counts the number of neighbor-ing particles in a shell with distance r to a particle. The function g110 (r) measures therelative orientation of two particles, where positive values represent parallel and negativevalues represent antiparallel average orientations. Finally, the function g101 (r) reflects theorientation of one particle with respect to the connecting axis pointing from the neighbor-ing particles to itself. Specifically, positive values of g101 (r) indicate that the hydrophilicside of a particle points away from the second particle. We also consider two higher-ordercoefficients, which are

g202 (r) ≡⟨

V

8πr2N2

∑

i

∑

j 6=i

δ (r − rij)(

3 (ui · rij)2 − 1)

⟩

(6.8)

and

g220 (r) ≡⟨

V

8πr2N2

∑

i

∑

j 6=i

δ (r − rij)(

3 (ui · uj)2 − 1)

⟩

. (6.9)

We see, that both Eq. (6.8) and (6.9) are closely related to Eq. (6.7) and (6.6), respec-tively. One difference is, that Eq. (6.8) and (6.9) do not distinguish the sign of the scalarproduct. Therefore, they can give additional information. If, e.g., Eq. (6.6) is zero, par-ticle orientations can be aligned perpendicular or parallel and antiparallel for the sum tovanish. Now, Eq. (6.9) can solve this uncertainty and in an analogous way Eq. (6.8) helpsto determine uncertainties of Eq. (6.7). Higher-order coefficients gl1l2,l are not consideredhere, since they do not give a better insight into the main feature of the system, thatis, aggregate formation. In fact, we will see later that in our simulations Eq. (6.8) andEq. (6.9) are not particularly important. The calculation of such pair correlation functionsin molecular dynamics simulations can be found in Ref. [77]. The Dirac delta distributionsare taken into account by the distance check of two particles. The calculation of any scalarproduct is straight forward and the normalization is done by calculation of the number ofideal gas particles, which would be found in the same shell. The volume of such a shell isV shell = 4π∆r((k + 1)3 − k3)/3, where k is a step index of discretized values rk and ∆r isthe width of the shell.

46

6.2. Cluster related order parameters and distributions

For our density functional theory study we use slightly different correlation functions, asgiven in Ref. [28]

h (z) =

∫

d2ωα (z, θ) u (ω) ez (6.10)

=2π

∫ π

0dθsinθα (z, θ) cosθ =

(

1

tanh (|a (z)|) − 1

|a (z)|

)

a (z) · ez

|a (z)| , (6.11)

h2 (z) =1

2

∫

d2ωα (z, θ)(

3 (u (ω) ez)2 − 1)

(6.12)

=π

∫ π

0dθsinθα (z, θ)

(

3cos2θ − 1)

= 3

(

1

|a (z)|2− 1

|a (z)| · tanh (|a (z)|)

)

+ 1,

(6.13)

where α (z, θ) is described by Eq. (4.90). The function h(z) describes the “polarization”.It can take values between −1 and 1 corresponding to complete alignment in negative andpositive z-direction for the considered position z. The latter, second-rank order parameterh2(z) also describes the ordering along the z-axis, but it is independent of the direction.Within our DFT calculations, the function h2(z) is strongly coupled to h(z), that is, h(z) =0 automatically yields h2(z) = 0 and vice versa. Hence, we mainly consider the z-dependent“polarization” h(z). In our discussion of results we show one comparison of both functions inFig. 8.3. In fact, both h(z) and h2(z) are somewhat related to their equivalent [cf. Eqs. (6.7)and (6.8)] in the spatial direction z, where the second particle position would be the wallposition and with a different normalization.

6.2. Cluster related order parameters and distributions

In the spirit of studies of chain formation [82], we also consider the fraction of boundparticles versus the number of particles N in the system. In our bulk simulations aggregateswith more than three particles form the dominant structures. We therefore redefine thisfraction as

φ =

⟨

∑

s≥4 sNC (s)

N

⟩

, (6.14)

where NC (s) is the cluster size distribution, which denotes the number of aggregates con-sisting of s particles, where s is the cluster size. Another interesting property is the meancluster size, which we define as the arithmetic average over all cluster sizes sk,

n =

⟨

1

N totC

NtotC∑

k=1

sk

⟩

=

⟨

1

N totC

∞∑

s=1

sNC (s)

⟩

=N

〈N totC 〉 , (6.15)

where N totC =

∑∞s=1NC (s) is the number of clusters. The second part of Eq. (6.15) gives

the relation to the previously defined cluster size distribution NC (s), and in the last termwe use the conservation of the particle number.

47


In addition, we define an order parameter for the description of aggregates with more orless spherical shape

GCluster101 =

⟨

1

N totC

∑

j

∑

i ǫ cluster j

ui · rij

Nj

⟩

, (6.16)

where rij is the unit vector from the center of cluster j to particle i and Nj denotes thenumber of particles in cluster j. This quantity is one for perfectly spherical clusters.

In fact, we need to define and find clusters at first, before we can evaluate these orderparameters and correlation functions. We explain our cluster search algorithm and theunderlying clustering criterion in Appx. C.1.

6.3. Time-dependent correlation functions

Assuming an equilibrated bulk system, where the system is homogeneous in all spatial direc-tions. Here it is sufficient to determine autocorrelation functions for one direction in space.Of course freezing transitions can destroy this symmetry, but in our molecular dynamicssimulations of the bulk system we have found the following autocorrelation functions tobe equal (within numerical accuracy) for all spatial directions. As an arbitrary choice weconsider the x-components. We study the velocity autocorrelation function, which is

Cvxvx (t) =

⟨

∑Ni=1 vxi

(t) vxi(0)⟩

⟨

∑Ni=1 vxi

(0) vxi(0)⟩ . (6.17)

This function is normalized by its value at zero time difference t, which can be derivedfrom the equipartition theorem for the kinetic energy and the temperature of the system.Further, we take a look at the orientational correlation function

Cuxux (t) =

⟨

∑Ni=1 uxi

(t) uxi(0)⟩

⟨

∑Ni=1 uxi

(0) uxi(0)⟩ . (6.18)

This function is also normalized by its value for a vanishing time difference. On longtimescales one expects both autocorrelation functions to vanish for typical fluids. Theappearance of negative values is understood as a cage effect for the velocity autocorrelationfunction, in which particles are trapped due to interactions and collisions with neighbors.Negative values of the orientational autocorrelation function are also caused by particleinteractions, which enforce a reorientation.

Since aggregation is an important characteristic feature of our systems, it is interestingto study the bond autocorrelation function,

Cbb (t) =

⟨

∑Ni=1

∑Nj=1,j 6=i bij (t) bij (0)

⟩

⟨

∑Ni=1

∑Nj=1,j 6=i bij (0) bij (0)

⟩ , (6.19)

where bij (t) is 1, if particles i and j are bound to each other, and otherwise it has thevalue 0. In a fluid system one expects bonds to break and to reform within a characteristic

48

6.4. Mean-square displacement

0.1 0.2 0.3 0.4

T*

0

1

2

3

4

5

6

7

8

9

10

C* V

ρ∗ = 0.01, ∼ρ∗ = 0.1, ∼ρ∗ = 0.3, ∼ρ∗ = 0.5, ∼ρ∗ = 0.6, ∼ρ∗ = 0.8, ∼ρ∗ = 0.01ρ∗ = 0.1ρ∗ = 0.3ρ∗ = 0.5ρ∗ = 0.6ρ∗ = 0.8

Figure 6.1.: Comparison of the specific heat capacity derived by Eq. (6.24) (denoted by ∼)with the one given in Eq. (6.23) as a function of the reduced temperature forvarious reduced densities, a detailed discussion of the results follows in Sec. 7.2.

timescale and therefore Cbb (t) should go to 0 for large times. We note that this bondautocorrelation function assumes an infinitely large system, otherwise one may subtract〈bij (t)〉2 [83].

6.4. Mean-square displacement

An interesting function in the context of dynamical behavior of fluids is the mean-squaredisplacement (MSD) [80], which is a measure for the covered distance from the originalposition, averaged over all particles in the system. It is defined as follows

⟨

∆r2 (t)⟩

=

⟨

1

N

N∑

i=1

[ri (t) − ri (0)]2⟩

. (6.20)

On short timescales one assumes a quasi undisturbed motion for typical fluids. This isnamed the ballistic regime and

⟨

∆r2 (t)⟩

∝ t2. When particle collisions and interactionsbecome prominent, which occurs on longer timescales, the diffusive (“Einstein-like”) regimeis found. It is characterized by

⟨

∆r2 (t)⟩

∝ Dt, where D is the diffusion constant.

6.5. Specific heat capacity

The occurrence of peaks in the specific heat capacity is an indication for phase transitions or,in our case, for aggregation. Therefore the specific heat capacity is an interesting quantity.

49


The heat capacity CV can be derived as follows [61, 76]

CV = (∂ 〈E〉 /∂T )V,N = − ∂

∂T

∂

∂βlnZN (T, V )

=kBβ2

(

1

ZN (T, V )

∂2ZN (T, V )

∂2β− 1

Z2N (T, V )

(

∂ZN (T, V )

∂β

)2)

=

⟨

E2⟩

− 〈E〉2

kBT 2. (6.21)

It is clear, that the contribution of the kinetic energy to the heat capacity is constant in thecanonical ensemble (the temperature is fixed), therefore we define the specific heat capacityas follows

CV =1

N(∂ 〈Epot〉 /∂T )V,N =

(⟨

E2pot

⟩

− 〈Epot〉2)

NkBT 2. (6.22)

For our simulations we consider this quantity in reduced units, so C∗V = CVk

−1B . As

mentioned in Sec. 5.2, the Berendsen thermostat does not reproduce the canonical ensemble.Morishita [79] derived a modified expression for the specific heat capacity, which is

CB∗V =

1

kBN(∂ 〈Epot〉 /∂T )V,N =

(⟨

E2pot

⟩

− 〈Epot〉2)

N (kBT )2 −B, (6.23)

where B is a correction given through

B =2

nN

(⟨

E2pot

⟩

− 〈Epot〉2)

√

√

√

√

√

⟨

E2kin

⟩

− 〈Ekin〉2

⟨

E2pot

⟩

− 〈Epot〉2. (6.24)

Here n is the number of degrees of freedom. This modified formula shows large derivations, ifthe energies Ekin or Epot fluctuate strongly. The fluctuation of the kinetic energy is coupledto an external heat bath through the coupling parameter τT of the Berendsen thermostat.In our molecular dynamics simulation we approximate the specific heat capacity usingEq. (6.23), because we are not interested in the value itself, but rather use the specific heatcapacity as a function of temperature as an indicator for aggregation. In Fig. 6.1 we showthe differences in the values of the specific heat capacity for both Eqs. (6.22) and (6.23)as a function of T ∗ for various ρ∗. The differences occur only in the peak regions, wherethe corrected values of C∗

V are larger, but do not shift the maxima. Therefore, it is stillpossible to consider Eq. (6.23) for our needs, but caution is needed when by comparing themagnitudes of C∗

V to other studies, CB∗V may be more appropriate.

6.6. Pressure

We will use the pressure to search for indicators of a liquid-vapor transition. Here wewill briefly derive a useful expression for our studies [in our DFT study we can derive the

50

6.6. Pressure

pressure from its relation to the grand canonical potential, Eq. (2.24)]. The equation of thepressure can be derived from the time average of the virial function of Clausius [69],

〈V 〉 ≡⟨

V(

rN)⟩

= limτ→∞

1

τ

∫ τ

0dt

N∑

i=1

ri (t) · ri (t)mi = −3NkBT, (6.25)

where V(

rN)

contains contributions of the particle interaction potentials V int and the

external potentials V ext acting on the N particles. In the last step of Eq. (6.25) we integrateby parts and use the equipartition theorem for the kinetic energy. Taking a closer look atthe contribution of V ext to Eq. (6.25) through the external walls, it follows that

⟨

V ext⟩

= −P∫

dSr · n = −P∫

∇ · rdV = −3PV, (6.26)

where dS is a surface element at position r, n is the normal vector of this surface elementand P is the pressure. By inserting Eq. (6.26) in Eq. (6.25) we find

P =NkBT

V+

⟨

∑Ni=1 ri · Fi

⟩

3V, (6.27)

where Fi is the sum over all forces exerted by other particles on particle i and V is thevolume.

51


52

7Chapter 7.

The bulk system of amphiphilic Janusparticles

3

This chapter focuses on the bulk system, i.e., a system without confining boundaries. Inthis way we can study the role of the bare particle interactions and the resulting aggregationbehavior. First, we start with a short DFT study in Sec. 7.1 which is followed by a MDsimulation study focusing on aggregation in Sec. 7.2. Further, we incorporate a slightlydifferent model for the MD study in Sec. 7.2.6 and compare the results with the formerones.

7.1. Density functional study of the bulk

Our first attempts to study the bulk system is motivated by the work of Tarazona andcoworkers [28]. They performed a similar DFT study for the present model and foundmembranes for a density distribution, which depends only on the Cartesian coordinate zand for only very dilute systems. We employ a different representation for the hard spherecontribution to the functional [cf. Eq. (4.89)]. In fact, the differences in the hard spherefunctionals are most pronounced in the vicinity of inhomogeneous structures, therefore acomparison of the results with respect to membrane formation would be interesting. Thefollowing results are also featured in our publication [1].

At high reduced temperature T ∗ and packing fraction η . 0.45, Eqs. (4.90) and (4.92)just yield a homogeneous, isotropic phase characterized by a constant number density,ρ (z) = ρ, and a constant orientational distribution, α (z, u (ω)) = 1/ (4π). Upon loweringthe temperature, inhomogeneous, “membrane-like” solutions can appear. The position ofsuch solutions are arbitrary, but we can pick out one realization by enforcing the position of“the membrane”. We refer hereby to an isolated bilayer of Janus particles, which is infinitelyextended along the x and y directions, and where the two sublayers are oriented towardeach other in the energetically preferable configuration (with facing hydrophobic sides). Todetect the corresponding temperature range in the present study, we start the numericalminimization with an initial guess for ρ (z), composed of two Gaussians. At high tempera-tures, this initial profile quickly disappears during the iteration procedure. For low valuesof T ∗, on the other hand, the procedure indeed yields the formation of periodically repeatedbilayer structures, characterized by a sharp decay of the density outside the bilayers anda depletion zone between the membranes. A typical density profile and a correspondingpolarization [Eq. (6.12)] as a function of the position z are shown in Figs. 7.1a) and b),respectively. The polarization in Fig. 7.1 b) shows that the bilayers are composed of two


53

7. The bulk system of amphiphilic Janus particles

58 60 62 64 66 68 70z / σ

-1

-0.5

0

0.5

1

h(z)

58 60 62 64 66 68z / σ

0

0.1

0.2

ρ(z)

σ3

a) b)

Figure 7.1.: Typical a) density profile ρ(z) and b) order parameter functions h(z) at η = 0.3and T ∗ = 0.2, when the bulk system starts to form bilayers for λ∗ = 3. In fact,due to growing bilayers and numerical errors [see thin peaks at z ≈ 63.5σ andz ≈ 69.5σ in a)] the density profile is not in equilibrium.

layers with particles with their hydrophilic side pointing away from one another. Further,the particles in the layers are completely polarized, as indicated by values of −1 and 1 ofh(z).

By determining the temperatures separating the homogeneous isotropic region from thebilayers for various packing fractions, we are able to draw the separation temperature line,which is drawn as the solid line in Fig. 7.2. The precise location of this line has to beconsidered with some caution since the inhomogeneous low-temperature profile does notconverge to a truly stable result, see Fig. 7.1. Therefore, one should consider the line inFig. 7.2 rather as an upper limit to the temperature range which allows membrane forma-tion. We also note that our results differ quantitatively from those obtained in Ref. [28],where the transition temperatures were found to be higher (although the same mean-fieldapproximation was used for the anisotropic interactions). These deviations might be dueto the fact that the HS part of the density functional used in Ref. [28] was treated in theweighted-density approximation, rather than by FMT, as in the present study. More im-portantly, in their (bulk) study Tarazona and coworkers allowed for spherical (particularlymicellar) structures in addition to the planar structures considered here. However, thetransition temperatures separating the homogeneous high-temperature phase from a phasewith free membranes, on one hand, and micelles, on the other hand, turned out to be quitesimilar [28]. For reasons discussed in Sec. 4.6 we do not search for spherical structures here.However, judging from Tarazona’s results we believe that our line plotted in Fig. 7.2 givesat least an estimate for the temperature below which ordered structure (of whatever type)exists. We recapitulate and discuss this point later when comparing these results with ourMD simulation study (cf. 8.2).

Apart from membrane formation, a further interesting question is whether a gas-liquidtransition can occur in our system. Such a transition was found, e.g., in a recent simulationstudy of patchy Janus particles [25]. Interestingly, the condensation in that system is“unconventional” in the sense that it is accompanied by pronounced cluster formation. Since

54

7.1. Density functional study of the bulk

0.1 0.2 0.3 η

0

0.1

0.2

T*

membranes

homogeneous

Figure 7.2.: Bulk phase diagram in the plane spanned by the reduced temperature T ∗ =kBTσ/C [where C is the strength of the Janus interaction appearing in Eq.(2.4)] and the packing fraction η = (π/6)ρσ3. The solid line indicates theestimated transition temperatures (see main text) between the homogeneous,isotropic high-temperature phase, and an ordered low temperature phase char-acterized by membranes. The dashed line corresponds to the vapor-liquid co-existence curve resulting from a density functional calculation with a MMFTtreatment of the anisotropic interactions (yet under the assumption of a ho-mogeneous isotropic state). Reprinted with permission from our study [1].Copyright 2011, American Institute of Physics.

our model is quite different from the one-patch Janus particles considered in Ref. [25] (seethe discussion at the end of Sec. 3.1), it is not immediately obvious that such a clustercondensation should occur in our system at all. However, considering the homogeneousisotropic system we do not find such an transition.

Nevertheless, we note that the absence of conventional condensation (within the homo-geneous phase) in our results is actually a consequence of the mean-field treatment ofthe anisotropic interactions; indeed, the corresponding term in the density functional [seeEq. (B.25)] is zero in a fully disordered system and therefore the system reduces to a systemof hard spheres. Interestingly, if we focus on such a homogeneous, isotropic system and usethe MMFT for the anisotropic contribution to the functional (see Sec. 4.6.1), we do finda gas-liquid phase transition, indicating that the Janus interaction is effectively attractive.The corresponding phase coexistence points are indicated by the dashed line in Fig. 7.2.Here, we do not explore this aspect further, the main reason being that our DFT studyfocuses on surface influences (see Sec. 8.1), where the use of MMFT in the presence ofsurfaces implies a drastic increase of numerical effort.

We discuss the bulk phase in the framework of molecular dynamics simulations (MD) inthe following section. Using these simulations, one is independent of any approximationsof interparticle correlations, as well as from an ansatz for the one-particle density.

55


0.1 0.2 0.3 0.4 0.5T

*0

1

2

3

4

5

6

7

8

9

10

C* V

ρ∗ = 0.01ρ∗ = 0.1ρ∗ = 0.3ρ∗ = 0.5ρ∗ = 0.6ρ∗ = 0.8

0.2 0.4 0.6 0.8 1T

*0

0.2

0.4

0.6

0.8

1

φ

a) b)

Figure 7.3.: a) Specific heat capacity (with some typical error bars, where the errors faraway from the maxima are negligible), b) fraction of particles in a cluster ofsize greater than 3 [see Eq. (6.14)] as a function of the reduced temperaturefor various densities.

7.2. Molecular dynamics study of the bulk

In the following we present MD results for the fluid phase regime ρ∗ . 0.8 and focus onaggregation behavior. Here, we follow our publication [2].

7.2.1. Aggregation diagram

Unlike our DFT calculations we can easily observe aggregation in MD simulations. Infact, in DFT the description of spherical aggregates may be included in the functional.Alternatively, as done by Tarazona and coworkers [28], one can consider one aggregateand fix its position, where it is possible to limit the spatial directions considered for thestudy of spherical aggregates to a radial-dependent density. In our MD simulations wedetermine the aggregation diagram and compare it with Fig. 7.2. First, we have to identifythe reduced temperature T ∗

agg = kBTaggσ/C, below which significant aggregation occurs. Infact, there is a variety of possible order parameters and thermodynamic quantities, whichcan be evaluated in this context. As a first indicator we consider the specific heat capacityCV at constant volume V (see Sec. 6.5). Results for C∗

V = CVk−1B as a function of the

reduced temperature are plotted for different densities in Fig. 7.3a). For each curve weobserve a pronounced maximum indicating an abrupt change of the potential energy in anarrow temperature range (or, equivalently, pronounced energy fluctuations in this interval).We therefore interpret the temperature related to the maximum as a measure of T ∗

agg. Theresulting aggregation line, which is obtained by plotting T ∗

agg as a function of the reduceddensity, is presented in Fig. 7.4 (solid line with squares). The advantage of using CV as anaggregation criterion is that this quantity is independent of any cluster search algorithm.

As alternative criteria for aggregation we focus on two quantities, which are based ondetermining the clusters. We note that the values of those quantities strongly depend on

56


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

ρ∗0.12

0.15

0.18

0.21

0.24

0.27

0.3

0.33

0.36

0.39

0.42

T*

T*

agg (φ = 0.5)

T*

agg (CV

= max)

T*

agg (Nc has min)

Figure 7.4.: Aggregation line, calculated from the specific heat capacity CV (solid line withsquares), the fraction of bound particles, φ, (dashed line), and from the clustersize distribution (solid line with triangles upside down). The lines connectingthe data points are guides for the eye. The indicated error bars are stronglyinfluenced by our discretization of the temperature. Additionally the error barsof the aggregation line resulting from CV are influenced by the error bars inFig. 7.3a). Reprinted with permission from our study [2]. Copyright 2012,American Institute of Physics.

the cluster search algorithm, which is discussed in detail in Appx. C.1. As a brief summarywe note that our clustering criterion counts two particles as connected if their distance isless than 1σ. If the distance is between 1σ and 1.2σ, the two particles belong to the samecluster, provided that they are orientated appropriately.

The quantity φ gives us a second option to find the aggregation temperature line using acluster search algorithm. It is defined as the number of particles bonded in a cluster of sizes greater than three, divided by the total number of particles [see Eq. (6.14)]. Results forφ as a function of reduced temperature and for representative reduced densities are shownin Fig. 7.3b). Using the functions φ (T ∗), one may define an aggregation temperature by

applying the relation φ(

T ∗agg

)

= 0.5 (see Ref. [82] for a similar strategy). Comparing the

two aggregation lines derived from CV and φ, see Fig. 7.4, we find similar results only atlow densities ρ∗ ≤ 0.1. At higher densities, the criterion based on φ predicts significantlyhigher aggregation temperatures. This is a consequence of our cluster search algorithmwhich, by definition, becomes difficult to apply for higher densities. In fact, it is not easyto estimate the error caused by this shortcoming of the cluster search algorithm. Therefore,we limit the use of this algorithm to ρ∗ ≤ 0.5. Still, there are probably some randomclusters incorporated at higher temperatures for ρ∗ ≤ 0.5, where member particles may benot bonded. We also note that the assumption that a cluster is composed of at least fourparticles [see Eq. (6.14)], is of course an arbitrary choice. This choice takes into accountthat we expect predominantly spherical cluster shapes at the densities considered. Indeed,we find that the inclusion of smaller structures into our aggregation criterion yields an evenless pronounced change of φ at T ∗

agg [as compared to the data shown in Fig. 7.3b)].

57


0 3 6 9 12 15 18s

0

10

20

30

40

50

60

NC

0 3 6 9 12 15 18s

0

10

20

30

40

50

60T

* = 0.25

T* = 0.2

T* = 0.14

T* = 0.11

T* = 0.1

0 3 6 9 12 15 18s

T* = 0.25

T* = 0.2

T* = 0.14

T* = 0.11

T* = 0.1

a) b) c)

Figure 7.5.: Number of clusters NC as a function of cluster size s for a) ρ∗ = 0.1, b) ρ∗ = 0.3,c) ρ∗ = 0.5.

To gain a third definition of the aggregation temperature, we consider the full clustersize distribution NC (s). Results for this function are plotted in Fig. 7.5 for representativetemperatures and densities. In a system where the cluster sizes are distributed randomly,the function NC (s) decays exponentially [84]. Therefore, the appearance of a minimumand a subsequent maximum at a specific cluster size can be interpreted as a signal ofaggregation [84]. Here we define Tagg as the temperature, at which the first minimum inthe function NC (s) appears (for each density considered). The resulting aggregation lineis also included in Fig. 7.4. Contrary to the aggregation temperature defined through φ,the definition utilizing NC (s) does not involve the appearance of random clusters. Thisexplains the better agreement with the first aggregation line (based on CV), particularlyat high densities. On the other hand, this third criterion does depend on the clustersearch algorithm, as does the one based on φ. Henceforth, we will take the aggregationline determined via NC (s) as our measure for aggregation as it is used in other works [84]and clearly more selective than the criterion based on φ. The comparison of Fig. 7.2 andFig. 7.4 shows that at densities ρ∗ ≤ 0.2 the aggregation line determined through NC (s)lies at significant higher temperatures than the temperature associated with the upper limitof membrane formation in the DFT.

As discussed in Sec. 7.1, a DFT treatment in the spirit of the MMFT shows the pos-sibility of a first-order phase transition, here a condensation. Triggered by this finding itis interesting to study the possibility of a condensation of a low-temperature aggregatedsystem, i.e., the transition between two states characterized by different cluster size distri-butions. Indeed, such a transition has been predicted in the simulation studies of Sciortinoand coworkers [25, 26, 53], too. The most promising way to search for such a transition isto perform grand-canonical Monte Carlo simulations (such as in Refs. [25, 26, 53]). Herewe restrict ourselves to MD simulations (incorporating a thermostat). To get at least a hint

58


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

ρ∗0

5

10

15

20

P*

T*=1

T*=0.2

T*=0.14

T*=0.1

a)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

ρ∗-30

-25

-20

-15

-10

-5

0

5

<E

* tot>

T*=1

T*=0.2

T*=0.14

T*=0.1

b)

Figure 7.6.: a) Pressure P ∗ and b) total energy 〈E∗tot〉 as a function of the density ρ∗ for

various temperatures T ∗. The definitions of the appearing reduced units aregiven in Appx. A.4.

for the presence or absence of a condensation transition, we investigate the pressure andthe potential energy as a function of density at various low temperatures and we plot someof the results of both quantities in Fig. 7.6. None of these functions reveals an indicationof condensation, such as a plateau or a loop-like behavior. This is a contrast to the resultsreported in Refs. [25, 26, 53], where a minimum in the function Epot (ρ) at intermediatedensities was found. Thus, the present Janus model system seems to lack a condensationtransition. The derivation of our DFT and MD simulation study seems to be caused bythe MMFT approach, which did not take the formation of clusters triggered by the highlydirectional character of the interaction into account. In fact, suppression of condensationdue to aggregation is also found in other models such as dipolar hard spheres [85].

7.2.2. Cluster properties

We have already shown the cluster size distributions in Fig. 7.5, which we reconsider in thissection. Further, we analyze the specific aggregate sizes and there properties found belowT ∗

agg in more detail. At this point we begin our discussion with the mean cluster size ngiven as Eq. (6.15). Results for n are plotted in Fig. 7.7. Upon lowering the temperatureclose to the aggregation points one observes a pronounced increase of n, which is similarto the behavior of φ [see Fig. 7.3b)]. Interestingly, at low temperatures the system withρ∗ = 0.01 seems to be characterized by somewhat larger mean clusters than the one atρ∗ = 0.1. However, these (small) differences may be due to equilibration problems in themost dilute systems. An illustration of typical clusters observed at temperatures close tothe aggregation line is given in Fig. 7.8. For small densities (ρ∗ . 0.3) and temperaturesclose above T ∗

agg, the system is dominated by small clusters containing up to 5 particles; seeFig. 7.8a) for two typical candidates. On the other hand, at lower temperatures one observes(in the same density range) aggregates of sizes 5 − 10, with a mean cluster size n ≈ 8 − 9.This behavior is also reflected by a corresponding peak in the cluster size distribution [seeFig. 7.5a)]. Similar average cluster sizes have also been observed in other studies involvingJanus particles at low densities, such as in the simulation studies of Sciortino and coworkers

59


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

T*

0

2

4

6

8

10

12

14

n

ρ∗ = 0.01ρ∗ = 0.1ρ∗ = 0.3ρ∗ = 0.5

Figure 7.7.: Average cluster size as a function of reduced temperature for various densi-ties. Reprinted with permission from our study [2]. Copyright 2012, AmericanInstitute of Physics.

[25, 26, 53]. Moreover, the cluster shapes found here resemble those in the experimentalstudy of Jiang et al. [23].

An further increase of ρ∗ toward the range ρ∗ = 0.3−0.6 shows that the low-temperaturesystem consists predominantly of clusters of size 13 [cf. Fig. 7.5c)], where for ρ∗ = 0.3 acompetition of such clusters of size 13 and smaller ones can be seen [cf. the double peakfor T ∗ = 0.1 in Fig. 7.5b)]. Surprisingly, we find the case T ∗ = 0.11 to differ for all threedensities in Fig. 7.5. For all other temperatures the trend appears to be clear, a reductionof temperature increases the size of clusters, but going from T ∗ = 0.11 to T ∗ = 0.1 we seethe reverse, which is almost not noticeable for ρ∗ = 0.5 and well pronounced for ρ∗ = 0.3.This is indeed even more puzzling than the finding that the trend of T ∗ > 0.11 to T ∗ = 0.1still fits. At this point we go back to clusters of size 13. As illustrated in Fig. 7.8c) (seefirst member in the row), these clusters are composed of one central particle surroundedby twelve neighboring particles. A closer inspection, including an analysis of correlationfunctions (see Sec. 7.2.3) shows that these 13-particle clusters are, in fact, icosahedrons[86, 87]. These are close-packed structures with five-fold symmetry. One main differencebetween an icosahedron and other arrangements involving twelve neighbors around a centralparticle, such as an fcc-cluster, is that the icosahedron cannot be periodically continuedto form an infinite, three dimensional lattice. On the other hand, it is known that theicosahedron maximizes the number of nearest-neighbor bondings. We suspect that it isthis “perfection of structure”, which leads to the overwhelming appearance of 13-particleicosahedral clusters in the system, as reflected by the delta-like peak in Fig. 7.5c). Wenote that one main feature of the icosahedron, that is, the appearance of a single particleas the cluster center, seems to be a unique property of the specific Janus pair interactionmodel used here [see Eq. (3.7)]. In particular as, according to our model, side-by-sideconfigurations with parallel orientations are energetically “neutral”. This is in contrast toother Janus potentials proposed in the literature [25, 26, 53]. Moreover, as discussed inSec. 7.2.6, the incorporation of terms into the anisotropic pair interaction potential which

60


Figure 7.8.: Representative clusters observed at a) temperatures above the aggregation line(T ∗

agg < T ∗ < 1), b) T ∗ < T ∗agg and densities ρ∗ . 0.3, c) T ∗ < T ∗

agg and ρ∗ & 0.3.In this context, the aggregation temperature is defined via NC (s). Reprintedwith permission from our study [2]. Copyright 2012, American Institute ofPhysics.

favor parallel orientations will prevent the encapsulation of a single particle within such astructure.

The formation of icosahedrons raises the question, what will happen upon a further in-crease of the density. At this point we reconsider Fig. 7.8c) and the cluster types shown andsee that at even higher densities (e.g., ρ∗ = 0.8) the icosahedrons have (at least partially)“merged” into even larger structures. As an example, one observes aggregates with morethan 30 particles composed of two icosahedrons linked by some sort of bridge. On theother hand, we do not find any vesicles or double-layer (membrane-like) structures. Suchstructures have been reported to occur at high densities, e.g. in the DFT study presentedin Ref. [28]. However, these authors included higher-order terms into the pair potential,Eq. (3.7). The impact of such additional terms in the context of our simulations will bediscussed in Sec. 7.2.6 and in Sec. 9. Here we merely note that vesicle-like structures(composed of 40 − 50 particles) have also been seen in another simulation study [26]. How-ever, that study involves yet another Janus model, which describes patchy particles withsquare-well attraction.

Interestingly, the aggregate structures observed in our system at high temperatures ratherresemble those reported in a MD study [88] of the melting and freezing of LJ clusters.Indeed, similar to our findings, the authors of Ref. [88] observe icosahedral clusters, whichcan interpenetrate to form huge polyicosahedrons [89].

7.2.3. Microscopic structure

Bearing the cluster size discussion of the previous section in mind, we now focus on the inter-nal structure of such clusters. Here, we begin by calculating some two-particle correlationfunction projections, i.e. Eqs. (6.5), (6.6), (6.7), (6.8) and (6.9). Higher-order coefficientsgl1l2,l are not considered here, since they do not give a better insight into the main feature

61


1 2 3r/σ

0

5

10

15

20

25

30

35

g 000(r

)

T* = 1

T* = 0.2

T* = 0.1

0.5 1 1.5 2 2.5r/σ

0

1

2

3

4

5

6

7

8

fccicosahedronT

* = 0.1

b)a)

Figure 7.9.: Radial correlation function at a) ρ∗ = 0.1, b) ρ∗ = 0.5 and various tempera-tures. Reprinted with permission from our study [2]. Copyright 2012, AmericanInstitute of Physics.

0 1 2 3 4r/σ

-2

-1

0

1

2

3

4

g 110(r

)

0 1 2 3 4 5r/σ

0

5

10

15

g 101(r

)

T* = 1

T* = 0.25

T* = 0.1

b)a)

Figure 7.10.: The correlation functions a) g110 (r) and b) g101 (r) at ρ∗ = 0.1 and varioustemperatures.

of the system, which is aggregate formation.

Numerical results for the functions g000 (r), g110 (r), g101 (r), g220 (r) and g202 (r) areplotted in Figs. 7.9 - 7.12, where we consider various temperatures around T ∗

agg. At hightemperatures (T ∗ & 1), the functions g110 (r), g101 (r), g220 (r) and g202 (r) are essentiallyzero everywhere (see Fig. 7.11), indicating that the local neighborhood of each particleis isotropic. At the same high temperatures, the function g000 (r) resembles that of asimple soft-sphere fluid. The aggregation process occurring upon lowering T ∗ (i.e., uponincrease of the Janus coupling) is already indicated at low densities [see Fig. 7.9a)] by thedevelopment of two sharp peaks at r ≈ 1σ and r ≈ 1.7σ. At the higher density ρ∗ = 0.5, thesecond peak has a significantly larger weight, reflecting the larger size of typical clusters and,consequently, the larger number of neighbors in the second shell. Indeed, as we have alreadydiscussed in Sec. 7.2.2, the dominant cluster size for ρ∗ = 0.5 and temperatures T ∗ ≈ 0.1 iss = 13. Moreover, in that section we have already suspected [based on the snapshot shown

62


1 2 3r/σ

0

1

2

3

g 110(r

)

1 2 3 4r/σ

0

1

2

3

g 101(r

)

T* = 1

T* = 0.25

T* = 0.1

b)a)

Figure 7.11.: The correlation functions a) g110 (r) and b) g101 (r) at ρ∗ = 0.5 and varioustemperatures. Reprinted with permission from our study [2]. Copyright 2012,American Institute of Physics.

0 1 2 3 4r/σ

-1.2

-0.8

-0.4

0

g 220(r

)

0 1 2 3 4 5r/σ

0

0.2

0.4

g 202(r

)

T* = 1

T* = 0.25

T* = 0.1

b)a)

Figure 7.12.: The correlation functions a) g220 (r) and b) g202 (r) at ρ∗ = 0.5 and varioustemperatures.

in Fig. 7.8c)] that the resulting cluster structure corresponds to an icosahedron, ratherthan to another, less densely packed, 13-particle structure such as an fcc-cluster. To clarifythis issue, we consider in Fig. 7.9b) the pair correlation function g000 (r12) at ρ∗ = 0.5 andT ∗ = 0.1 in more detail. The peak positions corresponding to an ideal icosahedron (whichcan be extracted from the lattice basis vectors of icosahedrons given, e.g., in Ref. [86]),as well as those corresponding to an ideal fcc-structure are also shown in Fig. 7.9b). Acharacteristic feature of the fcc-structure (but not the icosahedron) is the appearance of apeak at rfcc ≈

√2σ. Inspecting our numerical data in Fig. 7.9b), we find that our g000 (r12)

indeed lacks such a peak; instead, the second peak occurs at r ≈ 1.66σ, which is very closeto the value rico ≈ 1.7σ expected in an ideal icosahedron. In this sense, our microscopicdata fully confirm the icosahedral structure indicated by the snapshots.

We now turn back to the orientational correlation functions g110 (r) and g101 (r). FromFigs. 7.10 and 7.11 we see that both functions develop pronounced oscillations upon loweringT ∗. At intermediate temperatures (i.e., 0.2 ≤ T ∗ ≤ 0.25) the function g110 (r) is negative

63


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8T*

0

0.2

0.4

0.6

0.8

1

G10

1Clu

ster

ρ∗ = 0.01ρ∗ = 0.1ρ∗ = 0.3ρ∗ = 0.5ρ∗ = 0.6ρ∗ = 0.8

Figure 7.13.: Order parameter for the formation of spherical clusters, GCluster101 , as a function

of reduced temperature for various densities. Reprinted with permission fromour study [2]. Copyright 2012, American Institute of Physics.

0 2 4 6 8r/σ

0

1

2

3

4

5

g 000C

lust

er(r

)

0 1 2 3 4 5 6r/σ

0

1

2

3

4

5T

* = 1

T* = 0.25

T* = 0.1

b)a)

Figure 7.14.: Correlation function between two cluster centers at a) ρ∗ = 0.1, b) ρ∗ = 0.5.

for most distances r. This reflects that the typical aggregates at these states are still sosmall (s . 5) that the average angle between the orientation vectors of two neighboringparticles is larger than 90 degrees (on the average). Correspondingly, the function g101 (r)is mainly positive, indicating that the two neighbors tend to “face away” from each other.Finally, at very low temperatures (T ∗ = 0.1), the high peaks in g110 (r) and g101 (r) indicatethat the orientational structure in the cluster is very “sharp”, as already suggested by thesnapshots in Fig. 7.8b) and 7.8c). The majority of the nearest neighbor particles arealigned somewhat parallel, indicated by g110 (r) > 0. They point away from each other, asdetermined by g101 (r) > 0. For the “next-nearest” neighbors we find g101 (r) > 0 again, butg110 (r) < 0, i.e., particles have antiparallel orientation and face away from one another.

We consider the higher-order orientational correlation functions g220 (r) and g202 (r) onthe example of ρ∗ = 0.5, as plotted in Fig. 7.12, because we find mainly clusters of size 13 inthis region at low temperature as shown in Fig. 7.5c). We find the intermediate temperature

64


case T ∗ = 0.25 to resemble the general trend of the correlation functions at T ∗ = 0.1 forthe nearest neighbors and to vanish for larger distances, therefore we only discuss T ∗ = 0.1.The negative value of g220 (r) for the nearest neighbors indicates angles of the orientationalvectors of approximately 83 − 97, the comparison with g110 (r) fixes the average angle tomore than 83 and less than 90. The correlation function g202 (r) states, that one particlefaces away from its neighbors significantly for this distance. Considering the values at adistance r ≈ 1.66σ, where the second peak of Fig. 7.9b) is located, we find an average angleof more than 90 to less than 97, through g220 (r) and g110 (r). Summarizing, we find thesehigher-order correlation functions to fix the angle more precisely, which is most useful forangles around 90.

However, a detailed analysis of the low-temperature cluster structures on the basis ofthe correlation functions in general is somewhat difficult due to the interference of variouscluster sizes. Therefore, we additionally consider the order parameter GCluster

101 defined inEq. (6.16). This order parameter is one for a perfectly spherical aggregate [see, e.g., thetetrahedron in Fig. 7.8a)], in other words a micelle, where all hydrophilic sides are pointingoutwards from the cluster center. Results for GCluster

101 as a function of the reduced temper-ature at various reduced densities are presented in Fig. 7.13. At all densities considered,the order parameter approaches values close to one at sufficiently low values of T ∗, indi-cating again the energetic preference of micellar structures in our model. Not surprisingly,the most perfect spherical structures are found at low densities, where aggregates are wellseparated from one another.

Finally, we present results for the radial cluster-cluster correlation function, gCluster000 (r12)

(here, rij is the distance between the centers of mass) in Fig. 7.14. The statistical definitionof this function is analogous to Eq. (6.5). At both densities ρ∗ = 0.1 and 0.5 consideredin Fig. 7.14, we see a marked change of gCluster

000 (r12) upon lowering T ∗. In particular, thefirst peak shifts to larger and larger distances, reflecting the growing of the average clustersize (see Fig. 7.8). Moreover, the denser system is characterized by a significant range ofcorrelations extending over several cluster radii.

7.2.4. Dynamic properties

The real strength of MD simulations is the accessibility of time evolution, such as time cor-relation functions. Here, we focus on the dynamics resulting from the aggregation process.First, we derive the mean-square displacement (MSD) [see Eq. (6.20)] and the normalizedautocorrelation of velocities [see Eq. (6.17)] to investigate the translational dynamics. Tostudy the rotational dynamics, we have calculated the normalized autocorrelation functionof the single-particle orientation vectors [see Eq. (6.18)]. Due to the overall isotropy of thesystems, all Cartesian components of the orientation and velocity vectors exhibit identicaltime correlations on the average. It is therefore sufficient to focus on, say, the x-component.

Results for the time-dependent functions defined in Eq. (6.20) and (6.18) are plotted inFig. 7.15 for the densities ρ∗ = 0.1 and ρ∗ = 0.5, and two representative temperatures.The corresponding results for Eq. (6.17) are shown in Fig. 7.16 for the same parameters.Considering first the MSD, we see from Fig. 7.15a) that all systems reveal a short-timeballistic regime, as expected. However, at later times, “normal” translational dynamics,where the ballistic regime directly goes over into diffusion, only takes place at high temper-

65


0.01 1 100t*

10-4

10-2

100

102

104

<∆r

2 (t)>

0 5 10 15 20 25 30t*

-0.2

0

0.2

0.4

0.6

0.8

1

Cu xu x(t

)

T* = 1, ρ∗ = 0.1

T* = 0.1, ρ∗ = 0.1

T* = 1, ρ∗ = 0.5

T* = 0.1, ρ∗ = 0.5

2

1a) b)

Figure 7.15.: a) Translational mean-square displacement, b) autocorrelation function of thex-component of the orientation vector as a function of the reduced time t∗ (seeAppx. A.4) for various reduced temperatures and ρ∗ = 0.1, 0.5. Reprinted withpermission from our study [2]. Copyright 2012, American Institute of Physics.

atures. Below the aggregation line and, in particular, at larger densities [such as ρ∗ = 0.5in Fig. 7.15a)], the dynamics is characterized by pronounced sub-diffusive behavior as re-flected by the “plateau” in the MSD (i.e.,

⟨

∆r2 (t)⟩

∝ tα with α < 1). Clearly, the observedslowing-down of the translational motion is due to the aggregation and the resulting dy-namic trapping of the particles into clusters. We recall in this context that the dominantcluster type at the state point considered here is the icosahedron, [see Fig. 7.5c)] whichis regarded as a particularly stable aggregate already from the static (energetic) point ofview. Having this in mind, the results in Fig. 7.15a) show that the static stability of theicosahedron is accompanied by a long-lived cage effect.

At the same time, rotational motion slows down as well. This can be seen from Fig. 7.15b)containing data for the function Cuxux (t) [see Eq. (6.18)]. At high temperatures, thesteep descent of Cuxux (t) signals that time correlations between the individual orientationsquickly die out. This decay is followed by a short period of negative correlations, indicatingreorientation due to interactions with the (nearest) neighbors. At low temperatures, on theother hand, there are no negative correlations. Instead one observes, particularly at thehigher density ρ∗ = 0.5, that the initial orientation decays only very slowly. This behaviorreflects again the special properties, in this case the rigidity of orientations, of particlestrapped in the icosahedral aggregates.

At high temperatures we recognize for both densities considered that Cvxvx (t) decayscontinuously, where in the dilute system velocities are correlated for much longer due tofewer particle interactions. Aggregation decreases the time correlation of the velocitiesdrastically and additionally we find negative values for Cvxvx (t), most pronounced for shorttimescales. This is again a sign for a cage effect. Further, we find pronounced negativevalues for the case of ρ∗ = 0.5 at larger times.

Having studied structural properties of aggregation as well as their influence on dynamics,it is interesting to calculate the normalized bond correlation function, which is defined in

66


0

0.5

1

Cv xv x(t

)

T*=1

T*=0.1

0 1 2 3 4 5 6

t*

0

0.5

1

Cv xv x(t

)

T*=1

T*=0.1

a)

b)

Figure 7.16.: Autocorrelation function of the x-component of the velocity as a function oft∗ for various reduced temperatures and a) ρ∗ = 0.1, b) ρ∗ = 0.5.

Eq. (6.19). We show the corresponding results as a function of time for two different reducedtemperatures below the aggregation line at two different ρ∗ in Fig. 7.17. For all statesconsidered, the bond autocorrelation function displays stretched exponential behavior, i.e.,Cbb(t) = a0exp (−(t/a1)a2) where a2 < 1. At the (higher) temperature T ∗ = 0.17 therelaxation times a1 for both densities are rather short (with a1 being smallest at ρ∗ = 0.5),yielding a relatively fast decay of the autocorrelation function. However, at the lowertemperature T ∗ = 0.14 the relaxation times are significantly larger, reflecting the existenceof long-lived bonds. On the other hand, the stretching exponents a2 are larger than at T ∗ =0.17. Altogether, the slowest decay of Cbb(t) is found at ρ∗ = 0.5, i.e., the density relatedto the formation of icosahedrons. We note that the appearance of stretched exponentialdecay in a bond autocorrelation function is consistent with what has been found in otheraggregating systems, such as colloidal gels [90, 91]. Here, stretched exponential decay(with a2 ≈ 0.54) occurs in the autocorrelation of ternary bonds. Further examples arecertain undercooled liquids [92, 83], particularly water, where the hydrogen-bond correlationfunction is characterized by values of a2 in the range 0.6 ≤ a2 ≤ 0.73 [83].

In addition, we determine the cluster size distribution as a function of the size s for varioussimulation times to get an idea how icosahedrons are formed at ρ∗ = 0.5 for T ∗ = 0.1. Theresults are plotted in Fig. 7.18. We find a formation of clusters of sizes 6 to 19 in thebeginning and at longer times we find clusters of sizes 11, 13 and 14 to dominate. Hence,a continuously recombination of larger and smaller clusters to clusters of intermediatesizes takes place for increasing simulation time which results in mainly icosahedrons atequilibrium. We stress that the equilibrium solution is on a large timescale. One must becautious with these cluster size distributions, which clearly depend on the initial conditions(cf. Appx. C.4) as long as equilibrium is not reached.

67


0.01 0.1 1 10 100

t*

0

0.2

0.4

0.6

0.8

1

Cbb

(t)

ρ* = 0.1, T

* = 0.17

ρ* = 0.1, T

* = 0.14

ρ* = 0.5, T

* = 0.17

ρ* = 0.5, T

* = 0.14

Figure 7.17.: Bond autocorrelation function as a function of the reduced time at twodensities (ρ∗ = 0.1 and ρ∗ = 0.5) and two temperatures (T ∗ = 0.17 andT ∗ = 0.14). The black solid lines are fits according to a stretched exponentialfunction, i.e., Cbb(t) = a0exp (−(t/a1)a2). The fit parameters from bottomto top are (a0 = 0.71, a1 = 3, a2 = 0.526), (a0 = 0.56, a1 = 32, a2 = 0.56),(a0 = 0.64, a1 = 306, a2 = 0.726), and (a0 = 0.49, a1 = 736, a2 = 0.82).

7.2.5. The influence of the interaction range

So far we have considered the influence of the density and of the reduced temperature, i.e.,the anisotropic coupling strength, on aggregation. We do not perform a rigorous study ofthe influence of the inverse range parameter λ∗, but we carry out some test calculationsin addition to our former results for λ∗ = 3. It should be stressed that in principle evena condensation transition can be caused by changing the inverse interaction length. Forexample the patchy-particle model of Ref. [25] shows a condensation transition, where a longattractive interaction length due to a square-well potential with range 0.5σ is considered.But as mentioned earlier MD simulations are not the best method to track such a transition.

Therefore, we focus on aggregation, where we choose the two cases ρ∗ = 0.1 and ρ∗ = 0.5for T ∗ = 0.1. In fact, as our DFT study at walls (cf. Sec. 8.1.1) shows, the differences inthe z-dependent density and polarization profiles between λ∗ = 3 and λ∗ = 4 are small(see Fig. 8.4) and a value for λ∗ which is too large would at some point probably destroythe capability of forming aggregates. Tarazona and coworkers [28] stressed that λ∗ < 1.5results in an unphysical structure formation of tetra-layers, which are composed of a bilayerand an additional layer at each side. This structure is stabilized, e.g., by the attractionof the first and third layer, which is not found in experiments. Considering all these facts,we choose λ∗ = 2. This causes a longer anisotropic interaction range and gives rise to theformation of larger aggregates due to the additional stability of particles in more distantneighborhood, as shown in the two screenshots in Fig. 7.19. We plot the correspondingcluster size distributions in Fig. 7.20 a) and the radial correlation functions in Fig. 7.20 b).In comparison to the results for λ∗ = 3 (cf. Fig. 7.9) the maxima in the radial correlationfunctions appear to be slightly shifted to larger distances, while the second and third

68


0 2 4 6 8 10 12 14 16 18 20 22 24 26s

0

10

20

30

NC

t*= 0.05

t*= 0.5

t*= 2.5

t*= 10

t*= 40

t*= 80

t*= 100

equilibrated

Figure 7.18.: Number of clusters NC as a function of cluster size s at ρ∗ = 0.5, T ∗ = 0.1 andfor various simulation times t∗.

maxima are much more pronounced. In fact, the cluster size distribution changes in bothcases (cf. Fig. 7.5). We find a coexistence of cluster sizes 9 − 16 including icosahedrons atρ∗ = 0.1, where for the shorter interaction length sizes 6 − 10 dominate. In the case ofρ∗ = 0.5 icosahedrons are absent, here the sizes are larger than 18 with the last pronouncedmaxima at 36. For both densities the Eq. (6.16) gives values about 0.94 showing that theclusters are still almost spherical. We find for both densities in Fig. 7.20 b) peaks at thepeak positions of the icosahedrons as plotted in Fig. 7.9 b). This as well as the snapshotsshow, that the icosahedral structure is still dominant, i.e., for ρ∗ = 0.5 the larger aggregatesare formed by one icosahedron with additional particles attached outside or by some kindof fusion of such icosahedrons.

7.2.6. Results for the bulk system of modified Janus particles

Motivated by the question of the interplay of higher-order contributions to the pair potentialwith the interactions of the model studied so far, we consider the anisotropic potential,Eq. (3.9), without the contribution of Eq. (3.10). The reason of omitting Eq. (3.10) ison one hand to reduce the number of parameters in our model and on the other hand tobe closer to the study of Tarazona and coworkers [28]. Indeed, the parameter space isincreased significantly and we therefore limit our study to the case of q2 = 0.5, also theradial dependency as well as the sign of the interaction strength can be different in general.The most important difference between the potential defined in Eq. (3.9) (q1 = 0, q2 =0.5) and the “basic” model (q1 = q2 = 0) is that the extended potential favors parallelalignment, and thus can stabilize parallel side-by-side configurations. The chosen valueof q2 maintains the amphiphilic character of the original model. Three typical snapshotsobtained at low temperatures are shown in Fig. 7.21. Clearly, micellar structures formin the dilute system as shown in Fig. 7.21a). However, those micelles possess differentcluster size distributions and a more complex reduced temperature dependency than in the

69


Figure 7.19.: Snapshots from simulations with λ∗ = 2, a) ρ∗ = 0.1, T ∗ = 0.1, b) ρ∗ =0.5, T ∗ = 0.1.

“basic” model, see discussion below. The dense system is shown in Fig. 7.21b) and c). Instrong contrast to the basic model, we observe that lamellar (membrane-like) and worm-likestructures appear upon lowering the temperature.

We start our discussion by considering the specific heat capacity, which is plotted inFig. 7.22 as a function of the reduced temperature T ∗ and for different reduced densitiesρ∗. We compare this result with former results of the “basic” model [see Fig. 7.3a)]. First,we observe a dramatic shift of the maxima toward lower T ∗. This is probably a frustrationeffect induced by the competition of both anisotropic pair interactions. We know fromthe “basic” model that at intermediate T ∗ smaller aggregates form, but the additional pairinteraction counteracts the formation of particle pairs with hydrophilic sides pointing awayfrom one another. The effective attraction of such an alignment is indeed reduced by theadditional pair interaction. Clearly, different values of q2 cause different results, largervalues than 1 completely destroy the amphiphilic Janus character. On the other handnegative values favor the amphiphilic character as well as parallel in line alignments andthus penalize parallel side-by-side orientations. At this point we conclude that most valuesof q2 would simply destroy our model.

Looking at Fig. 7.22 we find hints for more than one maximum for most ρ∗. We find twosmall “peaks” in the dilute system, ρ∗ = 0.01. For ρ∗ = 0.1 only one “peak” is found. Atρ∗ = 0.3 and 0.5 the system shows more pronounced maxima than in the case of the “basic”model. Furthermore, the second “peaks” at lower T ∗ are indicated by the rising value ofC∗

V . With regard to our snapshots for the dense system (ρ∗ = 0.5) the first “peak” is likelyto relate to the appearance of bilayers and the second “peak” to the formation of worm-likestructures. Anyway, the occurrence of a second “peak” at lower densities is not explicablethat way. However, the number of T ∗ values is clearly not sufficient to track the completepeaks and their maxima.

We consider the cluster size distribution given in Fig. 7.23 as a function of the cluster sizes for two ρ∗ and various T ∗. At densities ρ∗ = 0.01 − 0.3 we find the formation of micelles,

70


0 5 10 15 20 25 30 35 40s

0

5

10

15

NC

0 1 2 3 4 5r/σ

0

2

4

6

8

10

12

14

g 000(r

)

ρ∗ = 0.1ρ∗ = 0.5

a) b)

Figure 7.20.: a) Number of clustersNC as a function of cluster size s and b) radial correlationfunction for two reduced densities and T ∗ = 0.1.

Figure 7.21.: Snapshots from simulations involving the extended potential, Eq. (3.9) (q1 =0, q2 = 0.5), a) ρ∗ = 0.01, T ∗ = 0.1, b) ρ∗ = 0.5, T ∗ = 0.05, c) ρ∗ = 0.5, T ∗ =0.1. We add a black line in b) to emphasize one worm-like structure.

but compared to the “basic” model we find structures of size s = 12 and additionallyat the two lower densities structures with s = 10 to be dominating at 0.1 ≤ T ∗ ≤ 0.07.A further reduction of T ∗ interestingly causes a reformation favoring smaller structures.As a first guess, we assume the influence of the additional interaction favoring parallelorientations to deform structures in this manner. Also, the occurrence of clusters of sizes = 12 contrary to those stable icosahedrons in the “basic” model shows the influenceof the additional pair interaction. In fact, considering snapshots [e.g. Fig. 7.21a)], thesestructures with s = 12 look almost identical to icosahedrons, but the center particle ismissing. Therefore, a further decrease of the temperature and with that an increase of theanisotropic coupling strength causes a deformation in less spherical clusters. This seemsto result in the breaking of structures into smaller ones, which are again spherical. Figure

71


0.05 0.1 0.15 0.2 0.25

T*

0

5

10

15

20

CV

*

ρ*=0.01 (add)

ρ*=0.1 (add)

ρ*=0.3 (add)

ρ*=0.5 (add)

ρ*=0.01

ρ*=0.1

ρ*=0.3

ρ*=0.5

Figure 7.22.: Specific heat capacity as a function of the reduced temperature for variousdensities. We plot the results of the extended model [indicated with (add)]and in addition, we show the old results of the original model [see Fig. 7.3a)].

0 4 8 12 16 20 24s

0

10

20

30

NC

4 8 12 16 20 24s

4 8 12 16 20 24s

10

20

30

T* = 0.2T* = 0.1T* = 0.07T* = 0.05

a) b) c)

Figure 7.23.: Number of clusters NC as a function of cluster size s for a) ρ∗ = 0.01, b)ρ∗ = 0.1, c) ρ∗ = 0.3.

7.24 shows the order parameter GCluster101 [see Eq. (7.13)] as a function of T ∗ for various

densities. GCluster101 gives a criterion to characterize the degree of sphericity. As expected,

the aggregates at T ∗ = 0.05 are more spherical than at higher T ∗ (the case ρ∗ = 0.01 beingpartly an exception).

In the dense case, ρ∗ = 0.5, only one or two aggregates form at low T ∗, including allparticles as shown in the snapshots, Fig. 7.21b) and c). These structures are connected totheir periodic image. For the discussion of the structure of the bilayers formed at T ∗ = 0.1we rotate the snapshot such that we can look on top of one of those layers (see Fig. 7.25)and omit the determination of an in-plane radial distribution function. The bilayer appears

72


0.05 0.1 0.15 0.2

T*

0

0.2

0.4

0.6

0.8

1

G10

1Clu

ster

ρ∗ = 0.01ρ∗ = 0.1ρ∗ = 0.3

Figure 7.24.: Order parameter for the formation of spherical clusters, GCluster101 , as a function

of reduced temperature for various densities.

Figure 7.25.: View from top on Fig. 7.21c), where both layers are emphasized by differentcolors.

to be composed of two almost quadratic layers, which are polarized such that the particlesin one layer are pointing away from particles in the other layer. The lattices of both layersare shifted diagonally relative to one another, i.e., a particle of the one layer has 4 neighborsof the other layer and 4 of its own.

73


74

8Chapter 8.

The influences of surfaces onamphiphilic Janus particles

4

Based on the study of the pure Janus particles in the bulk, we are now capable to discussthe influence of surfaces. Our aim is to explore the competition of surface-fluid, on theone side, and fluid-fluid interactions, on the other side. We begin our study within theframework of density functional theory in Sec. 8.1 and compare these results later with ourmolecular dynamics simulation in Sec. 8.2.

8.1. Density functional study of the influence of surfaces

We apply a pure mean-field DFT study in this section allowing only density profiles de-pending on z, the position normal to the surface, i.e., assuming translational invariance inthe other two spatial directions. The theoretical background is given in Sec. 4.6.

Reconsidering the bulk study, it is clear that the study of surface effects should beperformed at temperatures and packing fractions in the whole range above the solid line inFig. 7.2, i.e., where the bulk system is disordered (see Sec. 7.1). Otherwise the bulk systemwould be already structured by bilayers, which would cause numerical problems (cf. Sec. 7.1).The following discussion of our DFT results in the vicinity of surface influences orientatesclosely on our publication [1].

The remainder of this sections is structured as follows. At first, we consider a neutral(hard) planar wall Eq. (3.16) in Sec. 8.1.1. This discussion is then extended to hydropho-bic and hydrophilic surface fields Eq. (3.17) in Sec. 8.1.2 and finally we take a look onconfinement between two neutral planar walls in Sec. 8.1.3.

8.1.1. The neutral wall case

The hard (“neutral”) wall [see Eq. (3.16)] gives the most basic surface-fluid interactionpossible and is therefore a good start of our discussion. We find the differences between lowand high densities to change continuous, therefore we show only results for two densities,one for the dilute and one for a more dense system. In Fig. 8.1b), we present density profilesfor a dilute system characterized by η = 0.05 and different temperatures T ∗ ≥ T ∗

m, whereT ∗

m ≈ 0.079 is the temperature related to membrane formation in the bulk (see Fig. 7.2).The corresponding order parameter function h(z) [see Eq. (6.12)] is plotted in Fig. 8.1c).At T ∗ = 10, the density profile agrees within the numerical accuracy with that of a HS fluid,consistent with the negligible values of the order parameter. By decreasing the temperature,

4Selections from this chapter have been reprinted with permission from our study [1]. Copyright 2011,American Institute of Physics.

75

8. The influences of surfaces on amphiphilic Janus particles

0 1 2 3 4 5 6 7 8z / σ

-0.4

-0.2

0

0.2

h(z)

0 1 2 3 4 5 6 7z / σ

0

0.05

0.1

0.15

0.2

ρ(z)

σ3

T* = 0.079T* = 0.08T* = 0.09T* = 10

b) c)

Figure 8.1.: a) Sketch of the system at low temperature, b) density profiles ρ(z) and c)order parameter functions h(z) at η = 0.05 and various temperatures T ∗. Hereand in the subsequent figures, the hard wall is located at z = 0 and the rangeparameter in Eq. (3.8) is set to λ∗ = 3).

0 1 2 3 4 5 6 7z / σ

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

h(z)

0 1 2 3 4 5 6z / σ

0

0.5

1

1.5

2

2.5

3

ρ(z)

σ3

T* = 0.28T* = 0.3T* = 0.4T* = 10

b) c)

Figure 8.2.: a) Sketch of the system at low temperature, b) density profiles ρ(z) and c)order parameter functions h(z) at η = 0.3 and various temperatures T ∗.

configurations with neighboring hydrophobic sides of the particles become more and morefavorable, yielding orientational ordering of the particles as visible in Fig. 8.1c). At the sametime, the density profiles change from the typical behavior of a dilute HS system, wherethe density maximum occurs directly at the wall, toward a soft, loosely packed structurewhere most of the particles agglomerate somewhat away from the wall [see, e.g., data atT ∗ = 0.08 in Fig. 8.1b)]. Moreover, the distance between the first two layers indicatedby the two close maxima in ρ(z) is smaller than one particle diameter. Analyzing thecorresponding orientational order parameters we find that in the first layer (located atz ≈ 1σ), the particles are oriented such that their hydrophilic side points preferentiallyto the wall, whereas in the second layer (at z ≈ 1.8 − 2σ), they tend to orient in theopposite way. Thus, despite the low density considered, a significant degree of orientationalordering is already present. We sketch this finding in Fig. 8.1a). As expected, both thetranslational and orientational ordering becomes more pronounced upon an increase of thedensity. This is illustrated in Fig. 8.2 for the exemplary case η = 0.3. We sketch the low

76


0 1 2 3 4 5z / σ

0

0.05

0.1

0.15

0.2

h 2(z)

0 1 2 3 4 5 6 7 8z / σ

-0.4

-0.2

0

0.2

0.4

0.6

h(z)

T* = 0.3

T* = 0.4

T* = 0.6

T* = 10

a) b)

Figure 8.3.: a) Order parameter functions h(z) and b) order parameter functions h2(z) atη = 0.3 and various temperatures T ∗.

0 1 2 3 4 5 6 7 8z / σ

-0.4

-0.2

0

0.2

0.4

h(z)

0 1 2 3 4 5 6 7 8z / σ

0

0.5

1

1.5

2

2.5

3

ρ(z)

σ3

λ∗ = 2λ∗ = 3λ∗ = 4

Figure 8.4.: Density profiles ρ(z) and order parameter functions h(z) (inset) at η = 0.3 andT ∗ = 0.4 for various λ∗.

temperature behavior in Fig. 8.2a) and present the density and orientational profile in b)and c), respectively. At the lowest temperature considered, T ∗ = 0.28, the Janus particlesclose to the wall are almost completely arranged into a double layer, as reflected by thehigh values of the first, pronounced maxima in ρ(z) located at z = 0.5σ and z = 1.5σ.The double-layer formation is accompanied by a high degree of “polar” orientational order,as seen from the pronounced change of the function h(z) [see Fig. 8.2c)] from negativevalues at z ≈ 0.5σ (indicating that the hydrophilic side tend to point toward the wall) topositive values at z ≈ 1.5σ. Coming back to the density profile ρ(z) [see Fig. 8.2b)] wesee that at distances z beyond the bilayer at the wall, the third maximum appears only atz ≈ 2.7σ, indicating that the next layer is slightly shifted toward larger separations. Thisis a result of the repulsion between the hydrophilic sides of a particle in the second andthe one in the third layer. Indeed, the preferred orientation of the third layer becomesapparent from the negative value of h(z) in Fig. 8.2c). We plot h (z) again and compareit within h2 (z) [see Eq. (6.13)] in Fig. 8.3. As mentioned earlier in the definition of bothfunctions h (z) and h2 (z) are both zero at the same z positions as shown in Fig. 8.3a) andb), respectively. Further, h2 (z) does not differ between parallel and antiparallel orientation

77


0 1 2 3 4 5 6 7z / σ

0

5

10

15

20

25

ρ(z)

σ3

0 1 2 3 4 5 60

1

2 T* = 0.31

T* = 0.4

T* = 10

0 1 2 3 4 5 6 7 8z / σ

-1

-0.5

0

0.5

1

h(z)

b) c)

Figure 8.5.: a) Sketch at low temperature, b) density profiles ρ(z) and c) order parameterfunctions h(z) at η = 0.3 and various temperatures T ∗ at a hydrophilic wall[the parameters in Eq. (3.17) are set to C2 = C, λ∗

2 = 3]. The inset in part b)shows the density profile on an enlarged scale.

with respect to the z-axis. We can conclude that for our study, h2 (z) does not provideadditional informations, besides the case h (z) = h2 (z) = 0, i.e., particles are isotropicallydistributed at position z. We stress that all of these orientational ordering effects occur attemperatures where the corresponding bulk system is still homogeneous (see Fig. 7.2).

Within the framework of neutral walls, we have also considered the influence of the rangeparameter λ∗ of the Janus particle interaction [see Eq. (3.8)]. It turns out, however, thatthe latter has no crucial effect for the systems considered here. In Fig. 8.4 the densityprofiles and order parameter functions for the exemplary case T ∗ = 0.4 and η = 0.3 and forthree different λ∗ are shown. Upon increase of λ∗ (i.e., decrease of the interaction range)one merely observes a decrease of the density maxima close to the walls, and a damping ofthe oscillations at larger distances. Clearly, larger λ∗ give less interesting density profiles.Further, considering smaller values of λ∗ should be treated with caution, due to the fact,Tarazona and coworkers [28] found an unphysical formation of tetra-layers (composed ofa bilayer and an additional layer at each side attracted by the more distant layer of thebilayer) for λ∗ too small, and therefore they restricted λ∗ ≥ 1.5.

8.1.2. Influence of surface fields

Hydrophilic and hydrophobic surfaces enable us to enforce specific preferred orientations ofthe Janus particles close to a wall and extend our study beyond pure confinement effectsinduced by a neutral wall. These (locally) preferred orientations are modeled by Eq. (3.17).First, we consider a wall preferring the hydrophilic sides (such as silica [49]). This situationimplies that the surface field supports the orientation already found close to neutral walls(see Sec. 8.1.1). Exemplary density profiles and order parameters are shown in Fig. 8.5,where the bulk packing fraction η = 0.3, and the fluid-fluid and fluid-surface interactions(as well as the corresponding range parameters) are of the same magnitude. Due to thehydrophilic (and thus, supportive) character of the wall, the density in the contact layer,as well as the corresponding orientational order parameter, is even enhanced as compared

78


0 1 2 3 4 5 6 7z / σ

0

5

10

15

20

25

30

35

40

45

50

55

ρ(z)

σ3

0 1 2 30

0.2

0.4

0.6

0.8

1 T* = 0.28T* = 0.3T* = 0.4T* = 10

0 1 2 3 4 5 6 7 8z / σ

-0.5

0

0.5

1

h(z)

b) c)

Figure 8.6.: Same as Fig. 8.5, but for a hydrophobic wall characterized by C2 = −2C.

1 2 3z / σ

0

1

2

3

4

ρ(1) (z

)σ3

C2

*= 0

λ2∗ =2,C

2

*=0.5

λ2∗ =2,C

2

*=1

λ2∗ =2,C

2

*=2

λ2∗ =3,C

2

*=0.5

λ2∗ =3,C

2

*=1

λ2∗ =3,C

2

*=2

λ2∗ =4,C

2

*=0.5

λ2∗ =4,C

2

*=1

a)

1 2 3 4z/σ

-1

-0.5

0

0.5

1

h(z)

C2

*=0

λ2∗ =2,C

2

*=0.5

λ2∗ =2,C

2

*=1

λ2∗ =2,C

2

*=2

λ2∗ =3,C

2

*=0.5

λ2∗ =3,C

2

*=1

λ2∗ =3,C

2

*=2

λ2∗ =4,C

2

*=0.5

λ2∗ =4,C

2

*=1

b)

Figure 8.7.: a) Density profiles ρ(z) and b) order parameter functions h(z) at η = 0.275and T ∗ = 0.4 at a hydrophilic wall for various C∗

2 = C2/C and λ∗2.

to the corresponding system with neutral walls. This can be seen when we compare, e.g.,results for the temperature T ∗ = 0.31 in Figs 8.5b) and 8.5c) with corresponding ones atT ∗ = 0.3 in Fig. 8.2. Contrary the neutral-wall case, however, the second density maximumat the hydrophilic wall is much smaller than the first. We understand this behavior as aconsequence of the fact that in the second layer, the hydrophilic orientation preferred bythe wall competes with that dictated by the fluid-fluid interaction; the latter rather favorsthe hydrophobic side orienting toward the wall. The competition is also reflected by theasymmetric shape of the maximum in the function h(z). Finally, considering the third den-sity maximum and comparing with Fig. 8.2, we find that this (and the subsequent) peak(s)is higher and that the depletion effect (relative to the second layer) is less pronounced thanthat at a neutral wall. On the other hand, the density in the depletion zone is lower in thecase of the hydrophilic (than at a neutral) wall. We conclude that despite its short-rangecharacter [see Eq. (3.17)], the surface field is still effective even at fairly large separationsfrom the wall.

After the discussion of the hydrophilic wall we focus on the opposite case of a surfacepreferring the hydrophobic side of the particles. Here, completely different behavior is

79


1 2 3z / σ

0

1

2

3

4ρ(1

) (z)σ

3C

2

*= 0

λ2∗ =2,C

2

*=-0.5

λ2∗ =2,C

2

*=-1

λ2∗ =2,C

2

*=-2

λ2∗ =3,C

2

*=-1

λ2∗ =3,C

2

*=-2

λ2∗ =4,C

2

*=-0.5

λ2∗ =4,C

2

*=-1

λ2∗ =4,C

2

*=-2

b)

0.5 1 1.5 2 2.5 3 3.5 4z/σ

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

h(z)

C2

*=0

λ2∗ =2,C

2

*=-0.5

λ2∗ =2,C

2

*=-1

λ2∗ =2,C

2

*=-2

λ2∗ =3,C

2

*=-0.5

λ2∗ =3,C

2

*=-1

λ2∗ =3,C

2

*=-2

λ2∗ =4,C

2

*=-0.5

λ2∗ =4,C

2

*=-1

λ2∗ =4,C

2

*=-2

b)

Figure 8.8.: a) Density profiles ρ(z) and b) order parameter functions h(z) at η = 0.275and T ∗ = 0.4 at a hydrophobic wall for various C∗

2 = C2/C and λ∗2.

0.5 1 1.5 2z / σ

0

0.1

0.2

0.3

0.4

ρ(z)

σ3

C2

*=0

λ2∗ =2,C

2

*=0.05

λ2∗ =2,C

2

*=0.2

λ2∗ =3,C

2

*=0.05

λ2∗ =3,C

2

*=0.2

λ2∗ =4,C

2

*=0.05

λ2∗ =4,C

2

*=0.2

a)

0.5 1 1.5 2 2.5z/σ

-0.6

-0.4

-0.2

0

h(z)

C2

*=0

λ2∗ =2,C

2

*=0.05

λ2∗ =2,C

2

*=0.2

λ2∗ =3,C

2

*=0.05

λ2∗ =3,C

2

*=0.2

λ2∗ =4,C

2

*=0.05

λ2∗ =4,C

2

*=0.2

b)

Figure 8.9.: a) Density profiles ρ(z) and b) order parameter functions h(z) at η = 0.05 andT ∗ = 0.16 at a hydrophilic wall for various C∗

2 = C2/C and λ∗2.

observed. This situation is depicted in Fig. 8.6, where η = 0.3 and we choose, for thepurpose of illustration, a fluid-wall coupling parameter twice as large as that of the fluid-fluid interaction [i.e., C2 = −2C in Eq. (3.17)]. As a result of the dominant surface field,the bilayer formed at neutral walls (see Fig. 8.2) completely breaks down, and one ratherobserves the formation of a monolayer at sufficiently low temperatures. Moreover, thismonolayer is build by particles whose hydrophilic side points away from the wall, contraryto what has been observed before. The subsequent layer then has the reverse orientation andis shifted to slightly larger distances (z ≈ 1.65σ) as a result of the Janus repulsion betweenfirst and second layer. Only behind these first two layers, one observes the typical dampingof the oscillations in the density and orientation profiles, indicating that the surface-inducedstructures smear out away from the wall.

Clearly, all of these effects depend on the ratio between the fluid-wall and fluid-fluidcoupling parameter, C∗

2 = C2/C, and on the range parameter λ2. These influences on thedensity profile and polarization profile are shown for a slightly lower η = 0.275 in Fig. 8.7

80


0.5 1 1.5 2z/σ

0.1

0.2

ρ(z)

σ3

C2

*=0

λ2∗ =2,C

2

*=-0.05

λ2∗ =2,C

2

*=-0.2

λ2∗ =3,C

2

*=-0.05

λ2∗ =3,C

2

*=-0.2

λ2∗ =4,C

2

*=-0.05

λ2∗ =4,C

2

*=-0.2

a)

0.5 1 1.5 2 2.5z/σ

-0.2

0

0.2

0.4

h(z)

C2

*=0

λ2∗ =2,C

2

*=-0.05

λ2∗ =2,C

2

*=-0.2

λ2∗ =3,C

2

*=-0.05

λ2∗ =3,C

2

*=-0.2

λ2∗ =4,C

2

*=-0.05

λ2∗ =4,C

2

*=-0.2

b)

Figure 8.10.: a) Density profiles ρ(z) and b) order parameter functions h(z) at η = 0.05and T ∗ = 0.16 at a hydrophobic wall for various C∗

2 = C2/C and λ∗2.

(hydrophilic wall) and Fig. 8.8 (hydrophobic wall) for various C2/C and λ2 at the exemplarycase T ∗ = 0.4. Specifically, upon increase of |C2/C| one observes an increase of the extremaof both, ρ(z) and h(z), and an enhancement of the depletion areas. Similar effects emergeupon an increase of λ2, that is, a decrease of the range of the surface field. The influence ondensity and polarization of the fraction C∗

2 = C2/C seems to dominate. In 8.8b) we see theswitch of the preferred polarization close to a wall, when the hydrophobicity is increased.A weak hydrophobicity seems to reorientate only particles close to the wall, where a largerone increases the polarization of the more distant density peaks, too. So far we did notdiscuss the influence of surface fields on the dilute system. At this point we discuss the caseη = 0.05 for a hydrophilic and a hydrophobic wall for T ∗ = 0.16. The results are plottedin Figs. 8.9 and 8.10 analogous to those of the denser system. From these plots we find,that the overall dependency of the wall parameters does not change with the density. But,it is unlikely to compare the magnitudes for both densities, because the T ∗ regime wherethe anisotropic interaction becomes dominant compared to the pure hard-sphere repulsionis different.

8.1.3. Confinement between two neutral walls

In addition to the study of the influences of one wall it is interesting to analyze a confinedsystem. We consider a slit-pore geometry, i.e., our system is confined between two planarwalls. For simplicity, we focus on the case of neutral walls. For large wall separations Lz,one expects the structure at either surface to become decoupled from that at the othersurface, yielding bulk-like behavior [i.e., ρ(z) = ρbulk, h(z) = 0] in between the two walls.Indeed, we have explicitly checked that under such weakly confined conditions, the single-wall behavior discussed in Sec. 8.1.1 is recovered at each of the walls. For smaller Lz,pronounced confinement effects appear. Exemplary density and orientation profiles areplotted in Fig. 8.11b) and c), where Lz = 7σ (and η = 0.3). To illustrate this case, weadditionally sketch the system at low temperature close to each wall in Fig. 8.11a). At hightemperatures (i.e., close to the HS limit) the system is almost bulk-like in the center ofthe slit-pore. Decreasing the temperature we observe layer formation (which is typical for

81


0 1 2 3 4 5 6z / σ

0

1

2

3

ρ(z)

σ3

0 1 2 3 4 5 6 7z / σ

-0.5

0

0.5

h(z)

T* = 0.3

T* = 0.4

T* = 10

c)b)

Figure 8.11.: a) Sketch of the system at low temperature, b) density profiles ρ(z) and c)order parameter functions h(z) at η = 0.3 and various temperatures T ∗ inpresence of two walls located at z = 0σ and z = 7σ, yielding a separation ofLz = 7σ.

any confined system), accompanied by the development of orientational order particularlyclose to the walls. At the wall separation considered, the polarization at low temperaturesis asymmetric in the sense that particles in the left contact layer point antiparallel tothose in the right contact layer. This is consistent with the fact that the low-temperaturesystem consists of seven layers, a structure which allows for three full bilayer structures(composed of oppositely oriented particles as discussed in Sec. 8.1.1) plus one single layer.Two more examples are shown in Figs. 8.12 and 8.13 corresponding to the cases Lz = 3.5σand Lz = 3σ, respectively. In the first case, the high-temperature system is characterizedby three layers of particles, with the middle layer being rather thick [see Fig. 8.12a)]. Uponlowering the temperature, the anisotropic fluid-fluid interactions yields a splitting of themiddle peak, reflecting that particles in the middle layer tend to arrange in a “buckled”structure where neighboring particles are somewhat shifted to each other with respect tothe z-direction. The corresponding orientation profile in Fig. 8.12b) reveals that in thisbuckled middle layer, the particles arrange in an antiparallel way. Note that the resultingin-plane arrangement is not particularly unfavorable, since according to our model, side-by-side configurations are energetically “neutral” (see Fig. 3.1). By assuming this rathercomplex structure the system at Lz = 3.5σ overcomes frustration effects. Even strongerconfinement, as it is the case at Lz = 3σ (see Fig. 8.13), then yields three pronouncedlayers of particles, with the contact layers pointing in the opposite direction. However,contrary to the situation at Lz = 3.5σ, the order parameter directly at the position ofthe middle layer is zero. Only at slightly left or right of the center, one finds a preferredJanus-like orientation. We therefore regard this case as a frustrated system. The differentmicroscopic configurations appearing in the strongly confined systems in dependency ofLz give rise to pronounced oscillations of the normal pressure Pz [as calculated from thecontact theorem, see Eq. (4.88)]. The importance of this quantity stems from the fact thatit is experimentally accessible, e.g., by colloidal-probe atomic force microscopy [93]. Resultsfor Pz(Lz) at two temperatures are shown in Fig. 8.14. As it is typical for confined, densesystems of spherical particles, the oscillations have a period of about one particle diameter.

82

8.2. Molecular dynamics study in confinement

0 1 2 3z / σ

0

0.5

1

1.5

2

2.5

ρ(z)

σ3

T* = 0.3T* = 0.4T* = 10

0 1 2 3z / σ

-0.4

-0.2

0

0.2

0.4

h(z)

T* = 0.3T* = 0.4T* = 10

b)a)

Figure 8.12.: a) Density profiles ρ(z) and b) order parameter functions h(z) at η = 0.3and various temperatures T ∗. The walls are located at z = 0σ and z = 3.5σ.Reprinted with permission from our study [1]. Copyright 2011, AmericanInstitute of Physics.

Compared to the high-temperature situation, the anisotropic interactions not only stronglyenhance the amplitude of the oscillations, they also lead to slightly asymmetric peak shapesand to a shift of the oscillations. Indeed, at low temperatures, the maxima of Pz(Lz) occurat multiples of the particle diameter, consistent with the previously discussed frustrationeffects, e.g., at Lz = 3σ. On the other hand, the complex structure seen, e.g., at Lz = 3.5σ,corresponds to a minimum of the normal pressure curve.


Alternatively to the extended model defined in Eq. (3.9) and studied in Sec. 7.2.6 lamellarstructures can also be stabilized (within the original model) by planar solid surfaces, asshown before in our DFT study. This suggests that solid surfaces generally induce an in-teresting competition between spherical, micellar structures on the one hand, and planar,membrane-like structures on the other hand. In fact, our DFT calculations are not capableof describing micelles due to the assumption of translational invariance in two spatial direc-tions. Indeed, we already know from our molecular dynamics study of the bulk in Sec. 7.2that aggregation in mostly spherical clusters occurs and that even above the solid line inFig. 7.2 (for ρ∗ < 0.5). Motivated by this shortcoming of our DFT study we apply MDsimulations. In these simulations we consider the original model, Eqs. (3.12) and (3.7), andadd a confinement created by two soft walls, Eq. (3.20), with wall distance Lz = 10σ. Here,the wall distance measures the distance between the positions of infinitely large repulsionsof the walls. Contrary to hard walls the repulsion of walls described by Eq. (3.20) addition-ally reduces the accessible volume for the particles and raises the question of an effectivewall distance, which is clearly smaller than Lz = 10σ. Moreover, this effective wall distancedefines an average density of our system, which is used to compare our results with ourformer bulk MD simulations as well as with our DFT calculations. In a slit pore we define

83


0 1 2z / σ

0

0.5

1

1.5

2

2.5

3

ρ(z)

σ3

T* = 0.3T* = 0.4T* = 10

0 1 2 3z / σ

-0.4

-0.2

0

0.2

0.4

0.6

h(z)

T* = 0.3T* = 0.4T* = 10

a) b)

Figure 8.13.: Same as Fig. 8.12, but for walls located at z = 0σ and z = 3σ. Reprintedwith permission from our study [1]. Copyright 2011, American Institute ofPhysics.

2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8L

Z / σ

0.2

0.25

0.3

0.35

βPσ3

T* = 0.4T* = 10

Figure 8.14.: Normal pressure as a function of the wall separation for two temperatures T ∗.

Figure 8.15.: Snapshots from simulations involving soft walls with wall distance Lz = 10σ(effective 9σ) and a) ρ∗ = 0.111, T ∗ = 0.1, b) ρ∗ = 0.556, T ∗ = 0.1.

84


Figure 8.16.: Snapshots from simulations involving soft walls with wall distance Lz = 10σ(effective 9σ) and a) ρ∗ = 0.722, T ∗ = 0.3, b) ρ∗ = 0.722, T ∗ = 0.14, c)ρ∗ = 0.722, T ∗ = 0.1.

Figure 8.17.: View from top on the bilayers parallel to the walls with a slit pore length ofLz = 10σ (effective 9σ) for the two cases a) ρ∗ = 0.722, T ∗ = 0.14 and b)ρ∗ = 0.889, T ∗ = 0.1, where both layers are emphasized by different colors.

this average density for our MD simulation as

ρ∗ =σ3

Lz − σ

∫ Lz− σ2

σ2

dzρ (z) , (8.1)

where the z-dependent density ρ (z) has the dimension of the inverse of the volume and aneffective wall distance of Lz − σ is chosen. The reason behind our choice of the effectivewall distance is that at high temperatures the dilute system shows maxima of ρ (z) atpositions z ≈ 1 and z ≈ 9 [cf. Fig. 8.20a) in the following subsection], where the anisotropicpair interaction has minor influence. Another reason is that the wall-fluid interaction isbased on an average particle-particle interaction considering wall particles of diameter σ[cf. Eq. (3.20)], whereby we integrate over all wall particle centers.

In Figs. 8.15 and 8.16 five snapshots for different reduced densities and reduced temper-atures are shown. The first impression is that micelles are still the preferred aggregate forat least ρ∗ ≤ 0.556 at temperatures below the aggregation temperature, 0.18 ≤ T ∗

agg ≤ 0.21(cf. Fig. 8.15). At ρ∗ = 0.722 and T ∗ = 0.14 bilayers form close to the wall, which encapsu-

85


0.1 0.15 0.2 0.25

T*

5

10

15

CV

*

ρ*=0.0111

ρ*=0.111

ρ*=0.556

ρ*=0.722

ρ*=0.889

ρ*=0.01 (bulk)

ρ*=0.1 (bulk)

ρ*=0.5 (bulk)

ρ∗ =0.8 (bulk)

Figure 8.18.: Specific heat capacity as a function of the reduced temperature for variousdensities. The results indicated with (bulk) are taken from the original model[see Fig. 7.3a)].

0 5 10 15 20 25s

0

10

20

30

40

50

NC

0 20 40 60 80 100s

0

10

20

30

40

50

T* = 1T* = 0.2T* = 0.14T* = 0.1

a) b)

Figure 8.19.: Number of clusters NC as a function of cluster size s for a) ρ∗ = 0.111, b)ρ∗ = 0.556.

late micelles in the center between the walls [cf. Fig. 8.16b))]. At first sight these bilayersappear to be highly polarized with hydrophilic sides pointing outwards, where each layerappears to be hexagonal ordered with some distortions [cf. with snapshots in Figs. 8.16b)and 8.17a)]. For an even higher density ρ∗ = 0.888 at T ∗ = 0.1 we find the hexagonal orderof each layer to be intact as shown in Fig. 8.17b). A further finding of Fig. 8.17b) is, thatboth layers are shifted relative to one another. Considering ρ∗ = 0.722 for even lower T ∗

the snapshots show again micelles and no pronounced bilayers [cf. Fig. 8.16c)]. Indeed, ourorder parameter GCluster

101 [see Eq. (6.16)] indicates the existence of nonspherical structuresat intermediate temperatures and the predominance of micelles at lower temperatures byvalues of 0.38 and 0.78 at T ∗ = 0.14 and T ∗ = 0.1, respectively.

86


For a more detailed discussion we start by considering the specific heat capacity, whichis plotted in Fig. 8.18 as a function of the reduced temperature T ∗ and for various reduceddensities ρ∗. The comparison with the C∗

V resulting from the original model [see Fig. 7.3a)]shows less pronounced maxima. We note that our discretization in T ∗ is insufficient todetermine, whether the maxima of C∗

V are just shifted or really possess a smaller magnitude.A second note is that we consider somewhat different ρ∗ in our study of the confined system.Additionally, we find at ρ∗ = 0.556 a larger value of C∗

V for T ∗ = 0.14 compared to ρ∗ = 0.5and ρ∗ = 0.6 in the bulk, which indicates a surface induced effect, but again it is not clearif the width and height of the peak of C∗

V change both at the same time or only one ofthem.

To obtain some additional information we consider the cluster size distributions plottedin Fig. 8.19 at two different ρ∗ and various T ∗. At low densities the aggregates appear tobe of similar size compared to the bulk system [see Fig. 7.5a)]. Of course this similarity isnot surprising, because the influence of walls in the dilute system should be small as longas the wall distance is not too small. The cluster size distribution for the denser system atρ∗ = 0.556 and low temperatures is also similar to the bulk case shown in Fig. 7.5b). Wesee only a slight difference, which is the appearance of a small second maximum locatedat s = 26. However, the dominant clusters are still icosahedrons with s = 13. We donot show the results of cluster size distributions for ρ∗ = 0.722 and ρ∗ = 0.888, wherethe cluster search algorithm is partly not applicable due to the large density. Anyway, thesnapshot in Fig. 8.16b) shows bilayer formation close to both walls for low temperatures atρ∗ = 0.722 and micelles are located between both bilayers. Obviously, those micelles in thecenter are so close that our cluster search algorithm fails in separating these structures. Forρ∗ = 0.888 we find similar structure formation in snapshots. Contrary to ρ∗ = 0.888 we findfor ρ∗ = 0.722 bilayers only to be present at a small temperature range around T ∗ = 0.14.We state that the formation of aggregates is mostly influenced at densities ρ∗ ≥ 0.556 for awall distance Lz = 10σ (effective distance: 9σ).

8.2.1. Microscopic structure

We studied the influence of planar walls within our DFT approach and concentrated ondensity profiles and the order parameter function giving the polarization, both dependingon the distance from the wall. We found also the influence of two walls to decouple fromeach other at large wall separation Lz. The wall separation in this MD simulation studyis large, but not large enough to decouple the wall influences at low T ∗ and for high ρ∗.The DFT study showed the formation of bilayers close to the wall in dense systems and atlower densities kind of a loosely packed bilayer. Here, we will compare these results withour MD study. It should be stressed that the simulation is not grand canonical in contrastto the DFT study. Anyway, a comparison of the formation of structures is tempting andqualitatively possible, but for a quantitative comparison one should determine the chemicalpotential and with that the density within the slit pore in a grand canonical Monte Carlosimulation and use this result as a starting point for the MD simulation. Moreover, ourMD simulation incorporates a soft-wall as well as soft-sphere repulsion, where the DFTcalculation considered their hard equivalents. This again hinders the direct comparison ofdensities, the density in the MD simulation should be larger to involve the effective sphere

87


0 1 2 3 4 5 6 7 8 9 10z / σ

0

0.05

0.1

0.15

0.2

ρ(z)

σ3

T* = 0.1

T* = 0.2

T* = 1

T* = 0.1 (DFT)

T* = 0.2 (DFT)

T* = 1 (DFT)

a)

0 1 2 3 4 5 6 7 8 9 10z / σ

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

h(z)

T* = 0.1

T* = 0.2

T* = 0.1 (DFT)

T* = 0.2 (DFT)

b)

Figure 8.20.: a) Density profiles ρ(z) and b) order parameter functions h(z) at ρ∗ = 0.111and for various temperatures T ∗.

0 2 4 6 8 10z / σ

0

0.5

1

1.5

2

2.5

3

ρ(z)

σ3

T* = 0.1

T* = 0.25

T* = 1

T* = 0.3 (DFT)

T* = 1 (DFT)

a)

0 2 4 6 8 10z / σ

-0.5

0

0.5

h(z)

T* = 0.1

T* = 0.25

T* = 1

T* = 0.3 (DFT)

T* = 1 (DFT)

b)


0 1 2 3 4 5 6 7 8 9 10z / σ

0

1

2

3

4

5

ρ(z)

σ3

T* = 0.1

T* = 0.14

T* = 0.3

T* = 0.3 (DFT)

a)

0 1 2 3 4 5 6 7 8 9 10z / σ

-0.5

0

0.5

1

h(z)

T* = 0.1

T* = 0.14

T* = 0.3

T* = 0.3 (DFT)

b)


88


diameter, which is smaller. On the other hand, the soft-wall potential repulsion reachesinto the system and therefore reduces the effective wall distance.

In Figs. 8.20 - 8.22 we show the corresponding density profiles and polarizations [seeEq. (6.12)] at different ρ∗ for various T ∗. The following discussion considers only the leftwall, because the density and the polarization profiles at the right wall show the sameproperties. We compare our MD results in a slit pore of an effective width of 9σ withthose of a DFT calculation with a wall distance of 9σ. Clearly, the soft walls in our MDsimulation allow interpenetration and consequently confinement effects at high densitiesare less pronounced than in our DFT calculations [cf. 8.20a)]. For a density of ρ∗ = 0.111we have the advantage of comparing the whole temperature range of MD simulation andDFT study. For high temperatures we find a good agreement of MD and DFT calculations,where due to different density definitions the average density in the pore appears to beslightly larger in the MD simulation, while both polarization profiles are almost isotropic(not shown). Considering lower temperatures we find that the density profile resultingfrom the DFT changes significantly from the high temperature behavior for T ∗ = 0.1,while a large difference appears already for T ∗ = 0.2 in our MD simulation. The densitypeak positions in Fig. 8.20a) for our MD results are shifted to the center of the system by∆z ≈ 1σ. Here, aggregation takes place, thus the center of mass of a cluster can be foundin larger distance from a wall than one would expect a peak of unbound particles in thedensity profile. Contrary to this behavior we find for T ∗ = 0.1 that the DFT calculationshows a double peak in the density profile, the peaks are located at 1.2σ and 2σ. TheMD simulation on the other hand shows only a single maximum located around 2σ. Thereason for this significant difference between both methods is that the DFT does not allowspherical structures.

The polarization profiles of DFT and MD simulation in Fig. 8.20b) show similar trends, nosignificant shift is found. The MD system is stronger polarized. The results show, that theassumption of translational invariance in two spatial directions in the DFT underestimatesthe polarization and causes different shifts in the density profiles as well as a double peakstructure due to a missing description of micelles.

We turn to the higher density ρ∗ = 0.556. The density profiles of both methods appearsimilar at T ∗ = 1, but we find the peaks to be shifted to the wall in our MD simulation dueto a dense packing and a soft wall. Hence, it is not surprising that the density oscillationis more pronounced in the DFT [see Fig. 8.21a)]. This shift of peaks in the MD simulationis also noticeable in the polarization profile at high temperature, plotted in Fig. 8.21b).Moreover, the systems appears to be stronger polarized in the MD simulation close tothe wall, while the DFT shows a higher polarization more distant from the surface. Atlower temperature the density profiles do not show other differences between both methods.Contrary the MD simulation polarization profile shows a positive double peak behind thenegative polarization at the wall, which is less pronounced in the DFT result. The reasonis that “spherical” aggregates and not bilayers are found at low temperatures in the MDsimulation as discussed in the beginning. These polarization peaks indicate, that particlesclose to the wall are preferentially facing this wall with their hydrophilic side, while par-ticles in the second density peak (with respect to the wall) are mostly orientated in theopposite direction. For even lower T ∗ we do not show DFT results, since they would showa completely layered system as mentioned in our DFT study at walls. The MD simulation

89


shows in Fig. 8.21a) that density maxima and minima are much stronger pronounced atT ∗ = 0.1. Additional smaller peaks arise between the two peaks close to the wall and inlarger wall distance. In fact, this structure in the density profile is caused by packing ofprimarily clusters of size 13 [cf. Fig. 8.19b)]. Furthermore, the system appears to be highlypolarized by the packing of icosahedrons as shown in Fig. 8.21b).

As mentioned before the formation of bilayers close to the walls is found at ρ∗ = 0.722.Here, we focus only on lower temperatures, since we have already seen for ρ∗ = 0.556 thatthe high temperature behavior gives a relatively good agreement of MD simulation andDFT. For T ∗ = 0.3 we again find a shift of the peaks in the density and the polarizationprofiles in our MD simulation results (cf. Fig. 8.22) toward the wall. Hence, confinementeffects are much more pronounced for the DFT results, where the particles can not penetratethe walls. Moreover, the DFT results already show at this temperature the formation ofhighly polarized bilayers, where the density peak in the center reflects packing problemsdue to the confinement (cf. Sec. 8.1.3). By reduction of the temperature to T ∗ = 0.14in our MD simulation we find a much better agreement of DFT and MD results. At thispoint the MD result in Fig. 8.16b) indicates a bilayer close to each wall, with particles ineach layer of the bilayer pointing outward of the bilayer. Micelles are found around thecenter between the two walls. The polarization profile shows a large polarization not onlyfor the bilayers but also for the micelle region [cf. Fig. 8.22b)]. The density profile withinthe double layer and behind the layer goes to almost zero indicating that such structuresare close to crystallization [cf. Fig. 8.22a)]. The center between both walls shows a smalldensity in the density profile, while to the left and right each two peaks are found. Thepositions of these two peaks are determined by the center of masses of the micelles. Thesemicelles repel one another, which explains the lower density in the center of the slit pore.The repulsion between bilayers and micelles results in the depletion region between bilayerand micelles.

Thus the bilayer formation close to the walls predicted by the DFT can only be confirmedfor ρ∗ ≥ 0.722, where the layers appear to be almost crystallized as indicated by an averagedensity ρ (z) ≈ 0 close to the layers. We find a good agreement of density and polarizationprofiles for both methods close to the surfaces, when we compare T ∗ = 0.14 (MD) andT ∗ = 0.3 (DFT). Clearly, the micelles in the center are not predicted by the DFT and thuswe find large deviations between both methods for the profiles at the center of the slit pore.Interestingly a further reduction of the temperature to T ∗ = 0.1 for our MD simulationresults in different density and polarization profiles close to the walls for ρ∗ = 0.722. Atthis lowest temperature considered micelles dominate in the whole system and no bilayersare found. The formation of micelles by decreasing the temperature in a system alreadyincorporating bilayers shows an interesting competition between surface-induced packingeffects and anisotropic coupling strength for ρ∗ = 0.722. We do not know so far, whetherfor T ∗ < 0.1 micelles can be found in the dense system ρ∗ = 0.888 or if the bilayers closeto the wall persist.

Summarizing, we find that our DFT shows a good agreement with our MD simulationsat high temperatures and in the structural description of bilayers. So far, a completecomparison of density and temperature regions appears to be unlikely, as the wall separationis different. We note that a DFT result at ρ∗ = 0.611 for T ∗ = 0.3 (not shown) shows abetter agreement with the MD simulation for ρ∗ = 0.722 at the same temperature, i.e. there

90


are not yet sharp bilayers formed in the DFT. The next step is clearly to consider soft wallsfor the DFT calculation to reach a better comparability.

We stress that bilayers do not form in the bulk nor is the slit pore system completelydescribed by bilayers in our simulations.

91


92

9Chapter 9.

Conclusions and outlook5

This thesis is concerned with amphiphilic Janus particles. Their interaction is describedby the low-order contributions of an extension of a general pair-interaction potential oflinear molecules in spherical harmonics (see Sec. 3) as originally introduced by Tarazonaand coworkers [28]. This work is based on different methods, first, we use classical grandcanonical density functional theory (DFT) as published in our study [1] to describe thebulk system of amphiphilic Janus particles (see Sec. 7.1) and later on the competitionbetween surface-fluid and fluid-fluid interactions (see Sec. 8.1). The second method is themolecular dynamics (MD) simulation incorporating the Berendsen thermostat, where westudy the bulk system of these amphiphilic Janus particles in more detail with the mainfocus on aggregation (see Sec. 7.2). The majority of the MD simulation results are publishedin our study [2]. Additionally, we perform molecular dynamics simulations to determinethe influence of a model modification supporting parallel side-by-side configuration (seeSec. 7.2.6). Moreover, MD simulations are carried out to study the effects of confinedgeometry caused by two soft-walls (see Sec. 8.2), which yields a qualitative comparison toour former DFT results.

We discuss the results in two parts. First, we start with the analysis of the bulk sys-tem. Here, our mean-field DFT calculations show that in the bulk such membranes format low temperatures, where we concentrate on bilayer-like structures due to the assump-tion of invariance in two spatial directions. The reasons for our first DFT calculations inthe bulk are to determine the differences caused by our more sophisticated (as comparedto the Tarazona study [28]) treatment of the hard-sphere repulsion in the framework offundamental measure theory, and further, to estimate the reduced temperature, where thebulk system starts the formation of membranes. The later result is used to ensure that ourstudy of surface influences is carried out at parameters, where the bulk is still homogeneousand isotropic. Another open point is the presence of a vapor-liquid transition in the sys-tem, which is absent in a pure mean-field treatment of the anisotropic interaction. Withina modified mean-field approach where the pair correlation function is approximated by aBoltzmann factor, we find that there is, at least, a tendency for condensation. This is inqualitative agreement to what has been found in a recent simulation study [25] of a relatedJanus model system, where, however, the bulk condensation transition is accompanied bypronounced micellization. In the context of surface systems a bulk condensation transitioncould have important consequences, since it would enable a wetting transition (on top ofthe structures already observed). Therefore, it would be very interesting to improve thepresent density functional approach for the surface systems beyond the mean-field level inthe future. Apart from the MMFT approach, a further strategy would be to use a (second-


93

9. Conclusions and outlook

order) perturbation expansion of the excess free energy. In this context, the orientationalfree energy of the inhomogeneous surface system would be described by the anisotropic(direct) correlation functions of the homogeneous bulk system, which can, in turn, be cal-culated by integral equation techniques, such as hypernetted chain theory (see Refs. [94]and [95] for a similar approach in the context of inhomogeneous dipolar systems).

To overcome the computational limitations of our DFT study enforcing translational in-variance in two spatial directions, we perform MD simulations. As a first step we determinethe aggregation temperatures as a function of the density, using three different measuresbased on energy fluctuations or the cluster size distribution. It turns out that the vari-ous approaches yield consistent results only at low densities, whereas at higher densities,details of the cluster search algorithm strongly influence the (distribution-based) results.Irrespective of these details, our aggregation lines are located at somehow larger tempera-tures than those reported by Tarazona and coworkers in their original DFT study [28] ofthe model. Below the aggregation line, the system is dominated by specific cluster sizedistributions, which depend on the state considered. We analyze these clusters throughdirect visualization, by investigating various two-body correlation functions, and via a mi-cellar order parameter measuring the degree of spherical shape of the clusters. At low andintermediate densities, the system consists of small, weakly correlated, and nearly spheri-cal clusters containing 6 − 8 particles on the average (and with cluster shapes resemblingthose in experiments [23]). A characteristic feature of the present model at low tempera-tures and moderate densities (such as ρ∗ = 0.5) is the aggregation into icosahedrons, i.e.,close-packed clusters characterized by 13 particles and five-fold symmetry. We note thatsuch icosahedral clusters also play a prominent role in other systems such as metal melts.Due to aggregation of the metal atoms into icosahedrons, which are quasi-crystalline inthe sense that they cannot periodically be continued to form a regular lattice, the metalmelts can be cooled down to temperatures far below the melting point [96]. Upon furtherincrease of ρ∗ in our system, the icosahedrons become linked to one another and finallymerge into huge, elongated structures (such as icosahedrons linked by a “bridge”). On theother hand, we do not find any lamellar structures or vesicles (spherical, hollow structurescontaining curved bilayers), which have been seen in DFT [28] and simulation or integralequation studies [53, 25, 26]. In our opinion, the lack of lamellar structures is again a fea-ture of our model which does not favor side-by-side configurations with parallel alignments(which would stabilize bilayers), but rather micelles. Interestingly, we find no indication ofa condensation of these clusters (at least not in the temperature range considered). Accord-ing to our former modified mean-field DFT study of the isotropic, homogeneous bulk, theangle-averaged interaction between two of our particles is indeed sufficiently attractive tostabilize a condensation transition. However, this mean-field-like approach does not takeinto account the formation of clusters triggered by the highly directional character of theinteraction. In fact, suppression of condensation due to aggregation is also found in othermodels such as dipolar hard spheres [85]. Even, the absence of any hint of a condensationin our MD simulation is not a 100% reliable fact. Therefore, a study in the frameworkof grand canonical Monte Carlo simulation would be quite interesting, which is capableto determine such a condensation, if existent. We mentioned earlier that a condensationtransition was found in a recent simulated study of another Janus model [53, 25, 26]. Thismodel is essentially a one-patch model, where the (uniform) patches of two neighboring

94

particles are coupled through an attractive square-well potential with a significant longerinteraction length (0.5σ) compared to our study (long interaction range but with a fastdecay). In the present model the particles have essentially two “patches” with varying in-teracting strengths depending on the mutual orientation. Therefore, the absence of such acondensation may be explained by our shorter interaction length or by the additional repul-sion of the hydrophilic side of the particles in our model. The latter energetically penalizesa close packing of micelles next to each other. A future study of a possible condensationtransition using Monte Carlo simulation should address also the question of the interactionlength and thus clarify the difference of our model to the patchy-particles [25, 26, 53].

In addition to static properties, we study the mean-square-displacement (MSD) andvarious time-dependent correlation functions. In the MSD we find plateaus which arecaused by hindered diffusion due to long-lived aggregates. A cage effect is also found inthe velocity autocorrelation function below the aggregation temperature. Thus indicating astrong influence of aggregate formation on the dynamical behavior. The interesting questionof bond lifetimes is studied by the bond autocorrelation function, yielding long lifetimes ofaggregates at low temperatures.

A more detailed study of the time-dependent cluster formation as well as the clusterrecombination by breaking and fusion seems to be very interesting for future work. So farwe conclude from the time-resolved cluster size distribution that icosahedrons form on costof smaller and larger clusters, which form earlier. But still there are open questions.

Triggered by Tarazona and coworkers, as well as by the lack of membrane-like structuresin our findings, we incorporate some of the next order contributions to the anisotropicinteraction, i.e., we set q1 = 0 and q2 = 0.5 in Eq. (3.9). These terms were alreadysuggested in the original DFT study of the Janus model [28]. The limitation of takingnot all contributions to this order was used by Tarazona to be able to solve the functionalderivation with respect to α(r, u) analytically. Further, experiments [23] as well as thetheoretical study of the patchy particles [25, 26, 53] showed the occurrence of worm-like[23] and membrane-like [25, 26, 53] structures as well as vesicles [25, 26, 53] and largeraggregates [25, 26, 53]. In fact, our aim is to stabilize side-by-side orientations to allow forlarger structures, while maintaining the main features of amphiphilic Janus particles, whichare described by the “basic” model. Indeed, we find membranes and worm-like structuresin the dense system, where the dilute system forms micelles with different sizes than for the“basic” model. The additional pair interaction counteracts somewhat amphiphilic behaviorof our particles and therefore shifts the aggregation to lower reduced temperatures. A moredetailed future study of the interplay between incorporated higher order pair interactions,which are coupled by q1 and q2, may reveal a possible stabilization of vesicles.

The second part of this work focuses on effects of surfaces. Here we consider a funda-mental measure DFT calculation with a mean-field treatment of the contribution of theanisotropic pair interaction to the functional. One key finding of our study is that due tothe presence of a surface, significant translational and orientational ordering related to thebilayer formation occurs under conditions where the bulk system is still homogeneous andisotropic. Thus, the surfaces seem to strongly support the structure formation, even whenthis surface is just a neutral (hard) wall. Moreover, we show that the details of the inho-mogeneous structure at the wall can be “tuned” by varying the surface potential. Indeed,walls preferring the hydrophilic part tend to enhance the bilayer structure seen already

95


at neutral walls, whereas hydrophobic walls typically induce a competition between mono-and bilayer structures. Finally, we consider confinement effects emerging from the presenceof two planar (neutral) surfaces. It turns out that for specific wall separations, there arepronounced frustration effects stemming from the interplay between the fluid-surface po-tential, which prefers planar layering, and the fluid-fluid potential preferring bilayers withdepletion areas in between. This competition is also visible in an experimentally accessiblequantity, that is, the normal pressure as a function of wall separation.

To overcome the major drawback of our DFT study, namely the neglect of sphericalstructures, which are indeed important at temperatures below the aggregation line, weperform additional MD simulations in a slit pore with intermediate soft-wall distance Lz =10σ. While the comparison with the bulk MD simulation results shows similar aggregates inthe dilute system. We find membranes close to the walls in the densest case considered at lowand intermediate reduced temperatures. In the dilute system differences of the density andpolarization profiles perpendicular to the walls between the DFT and MD simulations canbe found. The finding of loosely packed bilayers in DFT seems to correspond to micelles inthe MD simulations, where the density and polarization profiles differ, because the internalstructure of the micelles is not involved in the DFT. Further, this internal structure causesa shift of the peaks of the density close to the walls according to the center of mass of themicelles. Nevertheless the profiles follow the same trend at high temperatures. Surprisingly,even for denser cases, the micelles are still dominant and show differences in the densityand polarization profiles of both methods. In the densest cases considered double layersclose to the walls are also found in the MD simulations, here the profiles are quite similar tothose of the DFT calculations close to the walls. Interestingly the MD simulation shows atρ∗ = 0.722 for T ∗ = 0.1 an absence of bilayers, which form at T ∗ = 0.14. Summarizing, ourpresent DFT results are more realistic under denser conditions (or at high temperatures) asexpected, while still micelles are not negligible in the MD simulations. Indeed, the strategyto include spherical self-assembled structures within our density functional approach toachieve a more realistic description of the system composed of Janus particles is generallyclear, as shown by Tarazona and coworkers [28] in their study of bulk systems. To achievea better comparability of MD simulation and DFT it would be interesting to incorporatesoft walls in the DFT to omit different effective slit pore sizes.

So far, most experiments on surface aggregation rather deal with true amphiphilic mole-cules. Nevertheless, one of the main effects reported in our study, namely, the fact thatthe surfaces strongly support the aggregation behavior relative to the bulk, is consistentwith the behavior of real surfactants (see, e.g., Ref. [97]). We also note that the changeof surface properties in our model can be experimentally realized by changing the degreeof hydrophobicity, e.g., by coating silicon wafers with polymer films of varying thickness[62]. Moreover, surface experiments involving dilute suspensions of elongated amphiphilicmolecules (rather than amphiphilic spheres) also report both planar and spherical structures[49, 62], suggesting that different structures could likewise occur for the spherical (Janus)case.

An interesting issue for future work concerns the impact of curvature of the substrateon the self-assembled structures [54, 98]. Indeed, from experiments [99] with amphiphilicmolecules it is well established that the curvature can be decisive for the type of structuresactually formed. Within such a study it would probably interesting to evaluate the interplay

96

of surface curvature and the choice of 0 ≤ q2 < 1 in Eq. (3.9) (q1 = 0). A comparison ofthe aggregates formed in the vicinity of such surfaces with those created by surfactantsof different chain length and head-to-tail ratios may be of interest, as our model can beinterpreted as a coarse-grained model of surfactants.

Given the recent progress in the synthesis of amphiphilic Janus nanoparticles of varyingcomplexity [29] we are confident that our theoretical predictions can, in principle, be testedby experiments. Anyway, the production of large quantities is still not accomplished [29],which probably will change in the future. The mentioned large variety of possible tech-nological applications for Janus particles in the introduction seems to be the beginning ofan exciting development in the field of such colloidal particles. Therefore, we assume thatJanus particles will be subject to theory and simulations in the future.

97


98

AAppendix A.

Appendix

A.1. Yukawa potential

The Yukawa potential can be motivated by a screened Coulomb potential. We consider anionic solution with electric neutrality, where the ions are not aligned periodically due to thesolution molecules weakening the Coulomb interactions and the thermodynamics. Ions aresurrounded by ions with opposite charge screening the center ion. Clearly, the higher thetemperature, the more uniform distributed are the ions. We assume a Maxwell-Boltzmannstatistic for the charge density

ρ (r) = ρ0 (r) e−qφ(r)(kBT )−1

, (A.1)

where φ (r) is the Coulomb potential and q is the charge. Considering a small perturbationin φ (r) caused by ρ1 (r) = qδ (r) it follows

φ (r) → φ (r) + δφ (r) . (A.2)

Assuming high T we can expand the exponential function in Eq. (A.1). Assuming, further-more, a uniform density (i.e. ρ0 (r) = ρ0), we obtain

ρ (r) ≈ ρ0

(

1 − qφ (r)

kBT− qδφ (r)

kBT

)

. (A.3)

The total change in the charge density is then approximatively given by

δρ (r) = qδ (r) − qρ0δφ (r)

kBT. (A.4)

At this point we use the Poisson’s equation,

∇2δφ (r) = −δρ (r)

ǫ0, (A.5)

with which follows that(

∇2 − λ2)

δφ (r) = −αδ (r) , (A.6)

where λ =√

qρ0/ (ǫ0kBT ) and α = q/ǫ0. Equation A.6 and incorporating the Fouriertransformation of φ (r) leads us to

(2π)− 32

∫

d3k(

k2 + λ2)

eik·rδφ (k) = αδ (r) . (A.7)

99

A. Appendix

Performing a further Fourier transformation of Eq. (A.7) yields

(2π)−3∫∫

d3rd3k(

k2 + λ2)

eir·(k−k′)δφ (k) = α (2π)− 32 . (A.8)

We now use the definition of the Dirac delta distribution, δ (k − k′) = (2π)−3 ∫ d3reir·(k−k′).It follows that

δφ(

k′) =α

(k′2 + λ2)(2π)− 3

2 . (A.9)

Performing a Fourier back-transformation and switching to spherical coordinates finallyyields

δφ (r) =α

2rπ2

∫ ∞

0dk

sin (kr)

(k2 − λ2)k (A.10)

=α′

re−λr, (A.11)

where α′ = α/ (4π).

A.2. Soft wall potential

The aim is to find the surface potential related to a thick unstructured planar soft wallwith a number density ρ. Therefore, we consider the soft-sphere potential Eq. (3.13) andintegrate over the wall volume V , i.e.,

φsoft wall = 4ǫwallρ

∫

Vd3r

σ12

|r − r′|12 . (A.12)

In the equation above we assumed the diameter σ of the wall particles to be equal to theone of the particles in our simulations. Further we limit ourselves to the one componentcase for the wall and for the fluid. The positions of the wall and the fluid particles arerepresented by r and r′, respectively. We apply cylindrical coordinates d3r = drrdxdφ withr in the range of 0 to ∞, x integrated from −∞ to xwall and φ in the interval of 0 to 2π.Further we substitute b ≡ b (x− x′) =

√r2 + x− x′ under the choice r′ = (x′, 0, 0)T , where

T is denoting the transposition. This yields

φsoft wall = −8ǫwallρ

∫ xwall

−∞dx

[

σ12

10b (x− x′)10

]b=∞

b=√

(x−x′)2

. (A.13)

In the next step we apply the upper and lower boundary for b, where the contribution for

the upper limit vanishes. Finally we substitute a =√

(x− x′)2 with x− x′ = −√a due to

the definition x′ > x and we arrive at the solution

φsoft wall (x) =4

45ρwallπǫwall

σ12

|x− xwall|9. (A.14)

100

A.3. Maxwell construction

0 0.1 0.2 0.3 0.4 0.5

ρ∗-0.01

0

0.01

0.02

0.03

0.04

0.05

P*

T*= 0.213

T*= 0.215

T*= 0.221

T*= 0.25

Figure A.1.: Pressure as a function of density for various temperatures, see Appx. A.4 forthe definition of the reduced units. The dashed lines should illustrate thecorresponding parallel lines for the Maxwell construction.

A.3. Maxwell construction

We follow the description of the Maxwell construction for the Van-der-Waals equationas given in Ref. [100]. Considering the pressure-volume isotherms a two-phase region isindicated by an occurrence of (∂P/∂V ) |T > 0 for temperatures below a critical temperature.The behavior (∂P/∂V ) |T > 0 is unphysical, because a volume reduction would cause apressure reduction. In this two-phase region the isotherm has to be replaced by a parallelline to the volume, V , axis. This is done by the Maxwell construction [100] in such a way,that the parallel line subdivides the maximum and the minimum such that both have thesame area. By this method we can construct the phase-coexistence curve. In fact, wecan also plot the pressure as a function of the density, an inverse measure of the volume,which is more useful for our DFT calculations. In Fig. A.1 we show the pressure curvesas a function of density for various temperatures. We have added the mentioned parallellines to emphasize what is meant with a subdivision into equal areas. The points of thecoexistence curve are then the outer intersection points of the parallel line with the pressurecurve. When both intersection points fall together, i.e., when the pressure curve has onlya saddle point, they define the critical point.

A.4. Reduced units

It is useful to define so-called reduced units in numerical calculations and simulations toomit discrete values and to offer a description which is valid for a wider branch of systems(e.g. with not specified diameter of the spheres, but ratios are fixed). The most important

101

A. Appendix

ones in our study are

t∗ = t√m−1ǫσ−1, (A.15)

T ∗ =kBTσ

C, (A.16)

ρ∗ = ρσ3, (A.17)

C∗V = CV k

−1B , (A.18)

λ∗ = λσ−1, (A.19)

P ∗ = σ3/ǫ, (A.20)

where t is the time, m is the mass, σ represents the sphere diameter, C is the anisotropiccoupling strength and λ is the corresponding (inverse) range (we use analogous definitionsfor C2 and λ2), ρ refers to the density, P is the pressure, T is the temperature, kB is theBoltzmann constant and CV is the specific heat capacity. Distances are used in units of σand the energy is treated in units of ǫ. The according reduced units of forces, torques andfurther are straightforward to derive and based on these given here.

A.5. Remark to the expansions of the pair potential and the paircorrelation function in rotational invariants

For our choice of correlation functions and the expansion of the interaction potential weneed the following spherical harmonics:

Y0,0 (ω) = (4π)− 12 , Y1,0 (ω) =

(

3

4π

) 12

cosθ, (A.21)

Y1,±1 (ω) = ∓(

3

8π

) 12

sin (θ) e±iφ, Y2,0 (ω) =

(

5

16π

) 12 (

3cos2 (θ) − 1)

,

(A.22)

Y2,±1 (ω) = ∓(

15

8π

) 12

sin (θ) cos (θ) e±iφ, Y2,±2 (ω) =

(

15

32π

) 12

sin2 (θ) e±i2φ. (A.23)

Further we find the occurring Clebsch-Gordon coefficients to be

C (0, 0, 0; 0, 0, 0) = C (2, 0, 2; 0, 0, 0) = C (1, 0, 1; 0, 0, 0) = 1, (A.24)

C (1, 1, 0; ±1,∓1, 0) = −C (1, 1, 0; 0, 0, 0) = 3− 12 , (A.25)

C (1, 0, 1; ±1,∓1, 0) = C (2, 0, 2; ±1,∓1, 0) = C (2, 0, 2; ±2,∓2, 0) = 0, (A.26)

C (2, 2, 0; 0, 0, 0) = C (2, 2, 0; ±2,∓2, 0) = −C (2, 2, 0; ±1,∓1, 0) = 5− 12 . (A.27)

In the final correlation functions we do omit the factor (4π)7/2 for (l1, l2, l) = (0, 0, 0),

(1, 0, 1) and (2, 0, 2), further, the factor − (4π)7/2 √3 for (1, 1, 0) and (4π)7/2 √

5 for (2, 2, 0)are neglected. The incorporation of u in spherical coordinates is done as follows, u1 · u2 =cos (θ1) cos (θ2)+sin (θ1) sin (θ2) cos (φ1 − φ2). The vector rij comes from the intermolecularframe with ω = (0, φ) such that rij is parallel to the z-axis.

102

BAppendix B.

Appendix

B.1. Proof for general introduction to DFT

B.1.1. Proof (1)

The proof [69] that f0 minimizes the functional Ω [f ] starts with Eq. (4.2) and appliance ofthe definition of f0 Eq. (4.1). This yields

Ω [f0] = Tr[

f0

(

HN − µN − β−1lnΞ −H + µN)]

= −β−1lnΞ = Ω. (B.1)

By reformulation of Eq. (4.1) we find a expression for HN − µN and can derive

Ω [f ] − Ω [f0] = β−1Tr (f lnf − f lnf0) . (B.2)

Considering the term

Tr (f lnf − f lnf0) = Tr

[

f0

(

f

f0ln

(

f

f0

)

− f

f0+ 1

)]

, (B.3)

we find it to be greater than zero, because of xlnx ≥ x− 1 for x > 0. In Eq. (B.3) we usethe normalization condition Eq. (4.2) and add Tr (f0) − Tr (f) = 0.

B.1.2. Proof (2)

The proof (based on Refs. [55] and [69]) is based on a contradiction of the initial assumption.This assumption is the existence of another potential V ′

ext, which results in the same densityρ0 (r,ω). We now consider the resulting Hamiltonian H ′

N , the resulting probability densityf ′

0 and the corresponding grand canonical potential Ω′. Ω′ and f ′0 refer to the original

temperature and chemical potential. Starting from Eq. (4.2) and applying Eq. (4.1) wefind that

Ω′ = Tr[

f ′0

(

H ′ −Nµ+ β−1lnf ′0

)]

≤ Tr[

f0

(

H ′ −Nµ+ β−1lnf0

)]

= Ω + Tr[

f0

(

V ′ − V)]

. (B.4)

Taking this into account, we use Eq. (4.9) yielding

Ω′ < Ω +

∫∫

d3rd3ωρ (r,ω)(

V ′ext (r,ω) − Vext (r,ω)

)

. (B.5)

By switching primed and unprimed expressions we also find that

Ω < Ω′ +

∫∫

d3rd3ωρ (r,ω)(

Vext (r,ω) − V ′ext (r,ω)

)

, (B.6)

103

B. Appendix

but then by adding both equations it follows that

Ω + Ω′ < Ω′ + Ω. (B.7)

Therefore, the assumption is wrong and there is only one Vext for a considered µ that givesrise to a given equilibrium density. Now we allow the µ′ in H ′ to be different from µ, thenthat u (r,ω) = µ− Vext (r,ω) is a unique functional of ρ0 (r,ω) follows. Further, it followsthan that f0 is a functional of ρ0 (r,ω) with f0 being a function of V .

B.1.3. Proof (3)

Following Ref. [69], we assume ρ (r,ω) and the corresponding f ′. We start as usual withEqs. (4.2) and (4.9) yielding

Ω[

f ′] =Tr[

f ′(

H −Nµ+ β−1f ′)]

=F in [ρ] +

∫∫

d3rd3ωρ (r,ω)Vext (r,ω) − µ

∫∫

d3rd3ωρ (r,ω)

=Ω [ρ] . (B.8)

With Ω [f ′] ≥ Ω [f0] and Ω [ρ0] = Ω [f0] = Ω it follows Ω [ρ0] ≤ Ω [ρ].

B.2. Numerical implementation of density functional theory6

In our study we carry out the minimization numerically, i.e., we solve Eq. (4.92) with aniterative algorithm [59] using a one-dimensional lattice with a total length of Lz = 60σ anda discretization of 512 − 1024 points per sphere diameter σ. Specifically, the density profilein step n is given by

ρn (z) = (1 − α) ρn−1 (z) + αeνρ0e−βφsurf,hard(z)−β 1

A

δFex[ρ]δρ(z)

∣

∣

ρn−1 , (B.9)

where νρ0 is defined in Eq. (B.43) and α is a mixing parameter 0 ≤ α ≤ 1 interpolatingbetween the old and new density profile. There are two constraints for the mixing parameterα. First, one must ensure in each step that n3 (z) < 1, because of the logarithmic term inEq. (B.48). Second, the convergence should be fast enough. In our calculations, we usethe following strategy (which partially follows that proposed in Ref. [59]): Within the firstiteration steps we calculate the grand potential for several values of the mixing parameter.The result is fitted (using a cubic fit) to find the value α0 where the functional becomesminimal. This value α0 is then used to define the new density distribution according toEq. (B.9). After few iterations with this procedure (each time updating α0), we keepthe value of α0 constant in the further iterations to minimize numerical noise. In thiswork α0 = 0.01αmax, whereby αmax controls the condition for n3 (z). As a criterion for

6Selections from this section have been reprinted with permission from our study [1]. Copyright 2011,American Institute of Physics.

104

B.2. Numerical implementation of density functional theory

convergence we monitor the difference between the old and new density distribution, thatis

∫

dz

∣

∣

∣

∣

eνρ0e−βφsurf,hard(z)−β 1

A

δFex[ρ]δρ(z)

∣

∣

ρn−1 − ρn−1 (z)

∣

∣

∣

∣

. (B.10)

The iteration is stopped when the integral in the equation above becomes less than 10−8.This method yields a better accuracy in terms of the boundary conditions, see Sec. 4.5.These boundary conditions provide a good estimation for the numerical accuracy. Withinour calculations, these boundaries conditions are fulfilled with very high accuracy (relativeerror less than 10−5) for all systems at neutral walls. In the presence of surface fields, thecontact theorem is fulfilled only with less accuracy (at worst 10−2 for the lowest temper-atures considered, we perform in addition calculations with a discretization of up to 8192and find a relative error at lowest temperature considered not exceeding 10−4), whereasthe bulk number density at large distance from the wall is still maintained. In the caseof a negligible amphiphilic interaction the adsorption at the wall is also a good controlparameter,

Γ =

∫ ∞

0dz (ρ (z) − ρbulk) = − ∂

∂βµ

(

∂φWBIIbulk

∂nbulk2

)

. (B.11)

For the numerical implementation we use algorithms from Ref. [101] and include the Gnuscientific library (see Appx. D).

A further important point is the initial density profile, we can increase the accuracyby starting with low densities, where the bulk number density as initial profile is a goodchoice. Then we reuse the final density profile as initial condition for a higher density. Alsoit is useful to decrease the reduced temperature stepwise and reinsert the results at highertemperature as initial profile. In the case of high anisotropic coupling strength this methodimproves the accuracy significantly.

To calculate the weighted densities, we rewrite Eq. (B.58) in Fourier space by using theconvolution theorem as emphasized in Ref. [59]. This yields

nβ (z) = FFT−1 (FFT (ρ (z)) · FFT (ωβ (z))) , (B.12)

where FFT stands for the fast Fourier transformation. The analogue of Eq. (B.12) holds forthe vectorial weighted densities nβ (z). In practice, the weight functions are multiplied withappropriate factors to allow for a polynomial interpolation within the integrals. Here weuse an order adaptive integration rule [102] with equivalently weighted internal nodes (11nodes) for every occurring integration. In the presents of a wall, one has to be careful withEq. (B.12) close to the wall. This is because the first interpolation factors are multipliedwith zero, here one has to correct the results. This is done for the calculation of theconvolutions in real space by using the order adaptive integration rule [102] wheneverpossible. Otherwise a polynomial of third order or, if even this was not possible, thetrapezoidal rule is applied.

In the occurring formulae are many expression which should be treated with caution.That is the necessity to calculate limiting values and to incorporate these with additional

105

B. Appendix

clauses. Otherwise the whole computation will fail or give rise to huge numerical problems,e.g. the expression

lim|a(z′+z)|→0

(

1

tanh |a (z′ + z)| − 1

|a (z′ + z)|

)

= 0 (B.13)

in Eq. (B.48) can course computational problems.

B.3. Functional derivation - how to do it

A more general introduction than the following short example can be found in Ref. [69].We assume a general functional G of a function f (x) such that

G [f ] =

∫

dnxg (f (x)) , (B.14)

where g is a function of f (x) and n represents the dimension of the vector x. The functionalderivation of G [f ] with respect to the function f is defined by [69]

δG [f ] =

∫

dnxδG [f ]

δf (x)δf (x) . (B.15)

Reconsidering Eq. (B.14) we variate the functional by increasing f (x) by δf (x), yielding

δG [f ] =

∫

dnxdg (f (x))

df (x)δf (x) . (B.16)

The comparison of this result with the former definition of the functional variation gives

δG [f ]

δf (x)=dg (f (x))

df (x). (B.17)

Further, we find that

δG [f ]

δf (x′)=

∫

dnxδG [f ]

δf (x)

δf (x)

δf (x′). (B.18)

Bringing all together we can derive

δG [f ]

δf (x′)=

∫

dnxdg (f (x))

df (x)

δf (x)

δf (x′)=

∫

dnxdg (f (x))

df (x)δ(

x − x′) =dg (f (x))

df (x)

∣

∣

∣

∣

∣

x=x′

,

(B.19)

in the third step of the equation above we use that we can reformulate f (x) as a functionalitself,

f (x) =

∫

dnx′f(

x′) δ(

x − x′) . (B.20)

106

B.4. Density functional theory calculations for Janus particles


B.4.1. The anisotropic contribution to the grand canonical functional

Within the mean-field approach we find the contribution of the anisotropic pair interactionto the functional to be

Ωan [ρ, α] =1

2

∫∫∫∫

d3r1d3r2d

2ω1d2ω2ρ (r1) ρ (r2)α (r1,ω1)

α (r2,ω2)φ1 (r12) (u (ω1) − u (ω2)) · r21. (B.21)

We then use the pair interaction Eq. (3.7) and that the unit vectors representing the orien-tation of a Janus particle, u (ωi), are subtracted from one another. Hence, we can write

Ωan [ρ, α] =1

2

∫∫

d3r1d3r2ρ (r1) ρ (r2)φ1 (r12) ·

·(

∫

d2ω2α (r2,ω2)

∫

d2ω1α (r1,ω1) u (ω1)

−∫

d2ω1α (r1,ω1)

∫

d2ω2α (r2,ω2) u (ω2)

)

· r21. (B.22)

After applying the normalization condition Eq. (4.13) we find that

Ωan [ρ, α] =1

2

∫∫

d3r1d3r2ρ (r1) ρ (r2)φ1 (r12) ·

·(∫

d2ω1α (r1,ω1) u (ω1) −∫

d2ω1α (r2,ω1) u (ω1)

)

· r21. (B.23)

Now we can use the fact that the exchange of particle one and two under the last integralcauses just a minus, yielding

Ωan [ρ, α] =1

2

∫∫∫

d3r1d3r2d

2ω1ρ (r1) ρ (r2)φ1 (r12) ·

· (α (r1,ω1) u (ω1) + α (r1,ω1) u (ω1)) · r21, (B.24)

where we obtain two similar expressions and this finally yields

Ωan [ρ, α] =

∫∫∫

d3r1d3r2d

2ω1ρ (r1) ρ (r2)φ1 (r12)α (r1,ω1) u (ω1) · r21. (B.25)

B.4.2. Functional minimization with respect to α

We remember the normalization condition Eq. (4.13), which gives us the functional variation

δΩ [ρ, α]

δα (r,ω)= λ (r) , (B.26)

107

B. Appendix

where λ (r) is the Lagrange multiplier enforcing Eq. (4.13). The functional derivation isthen

δΩ [ρ, α]

δα (r,ω)= β−1 δ

δα (r,ω)

∫∫

d3r1d2ω1ρ (r1)α (r1,ω1) ln (4πα (r1,ω1))

+δ

δα (r,ω)

∫∫

d3r1d2ω1ρ (r1)α (r1,ω1)φsurf (r1, u (ω1))

+δ

δα (r,ω)

∫∫∫

d3r1d3r2d

2ω1ρ (r1) ρ (r2)φ1 (r12) ·

· α (r1,ω1) u (ω1) · r21, (B.27)

which after renaming r = r1 yields

δΩ [ρ, α]

δα (r1,ω)=β−1ρ (r1) + β−1ρ (r1) ln (4πα (r1,ω)) + ρ (r1)φsurf (r1, u (ω))

+ ρ (r1)

∫

d3r2ρ (r2)φ1 (r12) u (ω) · r21. (B.28)

We reformulate the expression above to

α (r1,ω) = eβλ(r1)ρ(r1)

−ln(4π)−1−βφsurf(r1,u(ω))−β∫

d3r2ρ(r2)φ1(r12)u(ω)·r21 (B.29)

and put the result in the normalization condition, that is

∫

d2ωα (r1,ω) ≡ 1. (B.30)

This finally leads us to the result

α0 (r,ω) =exp (a (r) · u (ω))

∫

d2ω′exp (a (r) · u (ω′)), (B.31)

where

a (r1) = β

∫

d3r2ρ (r2)φ1 (r12) r12 − βφsurf (r1 · ez) . (B.32)

In the bulk case ρ (r) = ρ = const we find that

a = βρC

∫ ∞

σ

∫ 2π

0

∫ π

0drdφdθrsin (θ) e−λ(r−σ)

cos (φ) sin (θ)sin (φ) sin (θ)

cos (θ)

=

000

(B.33)

and it follows

α0 =1

4π. (B.34)

108


B.4.3. Functional minimization with respect to ρ

At first we derive the functional minimization for the bulk case with ρ (r) = ρ = const,

1

V

δΩ [ρ]

δρ

∣

∣

∣

ρ0

= 0, (B.35)

where V is the volume of the system. This result will help us to find the chemical potential,µ, which is related to a given bulk number density ρ. Later on we will control with the bulknumber density the chemical potential, µ = µ(ρ), in the case of incorporated walls. Thischoice for µ is taken, because the system will be in contact with the bulk and the anisotropicinteraction vanishes in the bulk within the mean-field approximation for a constant numberdensity. We now perform the functional derivation which reduces in the case of a constantnumber density to a derivation of a function, that is

1

V

δΩ [ρ]

δρ

∣

∣

∣

ρ0

=δ(

β−1ρ(

ln(

ρΛ3)

− 1)

− µρ+ β−1φWBII (nα))

δρ

∣

∣

∣

ρ0

(B.36)

= β−1ln (ρ0) + β−1ln(

Λ3)

− µ+ β−1 δφWBII (nα)

δρ

∣

∣

∣

ρ0

. (B.37)

The according weighted densities in the bulk have the equations

n3 =4πρR3

3, n2 = 4πρR2, n1 = ρR, (B.38)

n0 = ρ, n2 = n1 = 0. (B.39)

Now, we calculate at first the derivatives of Eqs. (4.73) and (4.74), which are

δψ2 (n3)

δρ=

4πR3

3n23

(

−2n3 − n23 − 2ln (1 − n3)

)

, (B.40)

δψ3 (n3)

δρ=

4πR3

3n33

(

−4n3 + 2n23 + 2n3

3 − 4 (1 − n3) ln (1 − n3))

. (B.41)

For the derivative of Eq. (4.72) we obtain

δφWBII (nα)

δρ= − ln (1 − n3) +

n3

1 − n3+

8πR3ρ

1 − n3

(

1 +ψ2 (n3)

3

)

+4πR3ρn3

(1 − n3)2

(

1 +ψ2 (n3)

3

)

+−2n3 − n2

3 − 2ln (1 − n3)

1 − n3

+

(

4πR2)3ρ2

8π (1 − n3)2

(

1 − ψ3 (n3)

3

)

+

(

4πR2)3ρ2n3

8π (1 − n3)3

(

1 − ψ3 (n3)

3

)

− −4n3 + 2n23 + 2n3

3 − 4 (1 − n3) ln (1 − n3)

2 (1 − n3)2 . (B.42)

Further, we define the quantity

νρ0 ≡ βµ− ln(

Λ3)

= ln (ρ0) +δφWBII (nα)

δρ

∣

∣

∣

ρ0

, (B.43)

109

B. Appendix

which we will use in our inhomogeneous functional derivation.

Putting α0 (r,ω) in the grand canonical functional then yields

Ω [ρ] =β−1∫

d3rρ (r) (ln (ρ (r)) − 1) − β−1νρ0

∫

d3rρ (r)

+ β−1∫

d3rφWBII (nα (r)) +

∫

d3rρ (r)φsurf,hard (r)

− β−1∫

d3rρ (r) ln

(

sinh (a (r))

|a (r)|

)

. (B.44)

In the following we assume a wall positioned on the z-axis with surface area A. Also,we assume translational invariance in x- and y-direction such that ρ (r) = ρ (z). Before wederive

1

A

δΩ [ρ]

δρ (z)

∣

∣

∣

ρ0(z)= 0, (B.45)

we calculate the expression

a (z) = 2πβezC

∫ σ

−σdz′ρ

(

z − z′) z′∫ ∞

2dbe−λσ(0.5b−1)

b− βφsurf (z)

= 2πβezCeλσΓ (0, λσ)

∫ σ

−σdz′ρ

(

z − z′) z′ − βφsurf (z) , (B.46)

which is derived by a convolution. The incomplete gamma function appearing in Eq. (B.46)is defined as

Γ (s, x) =

∫ ∞

xts−1e−tdt. (B.47)

The functional derivation then gives us

β

A

δΩ [ρ]

δρ (z)

∣

∣

∣

ρ0(z)= − νρ0 + ln (ρ) + βφsurf,hard (z)

+

∫ σ/2

−σ/2dz′(

ω2

(

z′)D2(

z′ + z)

+ ω2(

z′)D2(

z′ + z)

+ ω3(

z′)D3(

z′ + z)

)

− ln

(

sinh (|a (z)|)|a (z)|

)

− 4πCβeλσΓ (0, λσ)

∫ σ

−σdz′z′ρ

(

z′ + z)

·

·(

(

tanh(∣

∣a(

z′ + z)∣

∣

))−1 −∣

∣a(

z′ + z)∣

∣

−1) a (z′ + z) ez

|a (z′ + z)| , (B.48)

where the White Bear mark II and the last term are derived by correlations. The abbrevi-

110


ations used in the equation above are defined as follows

D2 (z) = − n2 (z)

2πR (1 − n3 (z))E1 (z) − 6n2 (z) n2 (z)

24π (1 − n3 (z))2E2 (z) , (B.49)

E1 (z) =

(

1 +

(

6n3 (z) − 3n23 (z) + 6 (1 − n3 (z)) ln (1 − n3 (z))

)

9n3 (z)

)

, (B.50)

E2 (z) =

1 −

(

6n3 (z) − 9n23 (z) + 6n3

3 (z) + 6 (1 − n3 (z))2 ln (1 − n3 (z)))

9n23 (z)

,

(B.51)

D2 (z) = − ln (1 − n3 (z))

4πR2+

n2 (z)

2πR (1 − n3 (z))E1 (z) +

3n22 (z) − 3n2

2 (z)

24π (1 − n3 (z))2E2 (z) ,

(B.52)

D3 (z) = − n2 (z)

4πR2 (1 − n3 (z))+

n22 (z) − n2

2 (z)

4πR (1 − n3 (z))2E1 (z)

+n3

2 (z) − 3n2 (z) n22 (z)

12π (1 − n3 (z))3 E2 (z)

− n22 (z) − n2

2 (z)

4πR (1 − n3 (z))

(

1

3+

2

3n3 (z)+

2ln (1 − n3)

3n23 (z)

)

+n3

2 (z) − 3n2 (z) n22 (z)

24π (1 − n3 (z))2

(

4

3n22 (z)

− 2

3− 4

3n23 (z)

(

1 − 1

n3 (z)

)

ln (1 − n3 (z))

− 2

3n3 (z)

)

. (B.53)

The assumption of translational invariance in two spatial directions influences the occurringweight functions, which then possess only z-dependency,

ω3 (z) = π(

R2 − z2)

, (B.54)

ω2 (z) = 2πR, (B.55)

ω2 (z) = 2πzez . (B.56)

The corresponding weighted densities are given as

n (z) =

∫ R

−Rdz′ρ

(

z − z′)ω(

z′) , (B.57)

n (z) =

∫ R

−Rdz′ρ

(

z − z′)ω(

z′) . (B.58)

111

B. Appendix

112

Figure C.1.: Flow diagram illustrating the cluster search algorithm. The algorithm startsby checking the distance rij between particles i and j, where j is a memberof cluster c. Reprinted with permission from our study [2]. Copyright 2012,American Institute of Physics.

CAppendix C.

Appendix

C.1. Cluster search algorithm and clustering criteria7

There are many different clustering criteria and cluster finding algorithms, reaching fromsimple distance dependent ones to more sophisticated ones like kmeans and kmeans++[103]. An example for a cluster criterion is the Davies-Bouldin index [104]. Probably thebest cluster search algorithm would take the interaction energy into account. But becauseof economic reasons we used the simple distance dependent version, introducing, however,some small modifications. Our algorithm is explained as a flow diagram in Fig. C.1. Wefirst evaluate the distance rij of two particles i and j. If they are closer than 1.2σ we doadditional calculations (see below); otherwise particle i will not join the cluster c, in whichparticle j is a member. On the other hand, if rij < 1.2σ particle i is included into thecluster if (at least) one at the following conditions is true i) the two particles overlap, ii)particle i is sufficiently close to the center of the cluster c, iii) particle i is pointing awayfrom the cluster center. We note that due to the dependency on the cluster center we runthe algorithm three times to unite clusters. The validity of this algorithm is checked bysnapshots and a comparison to a kmeans++ algorithm incorporating the Davies-Bouldin

7Selections from this section have been reprinted with permission from our study [2]. Copyright 2012,American Institute of Physics.

113

C. Appendix

index for different representative cases of our study.

In our molecular dynamics simulations including an additional anisotropic contribution tothe pair potential, Eq. (3.9) (q1 = 0), and in our simulations involving walls we have to dealwith bilayers and worm-like structures, too. Here, we cannot apply the clustering criterionabove. We therefore take only a distance dependent one, which puts each particle in acluster, if the distance to a particle in the cluster is smaller than the first minimum positionof the radial distribution function [see Eq. (6.5)]. Unfortunately this causes some failuresfor higher temperatures, where the first minimum in the radial distribution functional liesat larger distances, we did limit the distance therefore to a maximum of 1.8σ. Anyway, thismethod can also find random clusters, which are not connected by attraction.

C.2. Forces and torques for molecular dynamics simulations ofJanus particles

C.2.1. Basic Janus interaction

We consider the pair potential defined in Eq. (3.12) with a soft sphere repulsion, Eq. (3.13),instead of the hard-sphere contribution and derive the pair forces,

Fij = −Fji = −∇riφ (rij, ui, uj) (C.1)

=48ǫσ12

|rij|14 rij +Ce−λ(|rij |−σ)

|rij|2

((

2

|rij | + λ

)

(ui − uj) · rji

|rij | rij + ui − uj

)

. (C.2)

Based on the same pair interaction potential we derive the torques,

Tij = −ui × ∇uiφ (rij , ui, uj) (C.3)

= −Ce−λ(|rij |−σ)

|rij |2ui × rji, (C.4)

Tji =Ce−λ(|rij |−σ)

|rij |2uj × rji. (C.5)

C.2.2. Incorporation of higher-order pair interactions

The second additional higher-order contribution in Eq. (3.9), q2φ3 (r12)Q2 (r21, u1, u2), tothe pair potential causes additional pair forces. Assuming φ3 (r12) ≡ φ1 (r12) we derive

Faddij =

q2Ce−λ(|rij |−σ)

2 |rij|2

rij

(

−2 − 2

|rij | +(ui · rji)

2 + (uj · rji)2

|rij|2

(

9

|rij|+ 3λ

))

− 6

|rij | (ui (ui · rji) + uj (uj · rji))

. (C.6)

114

C.3. Cell lists and neighbor lists

Further, an additional contribution to the torques appears, that is

Taddij = −3q2Ce

−λ(|rij |−σ)

|rij|3(ui · rji) (ui × rji) , (C.7)

Taddji = −3q2Ce

−λ(|rij |−σ)

|rij|3(uj · rji) (uj × rji) . (C.8)

C.2.3. Forces at an unstructured planar soft wall

The incorporation of unstructured planar soft walls, Eq. (3.20), gives for each wall anadditional contribution to the, e.g., x-component of the force on particle i, that is

Fwallx,i =

36ρwallπǫwallσ12

45 |xi − xwall|11 (xi − xwall) . (C.9)

C.3. Cell lists and neighbor lists

We incorporate cell lists and neighbor lists of cells [76] to reduce the computation time.A cell is a box with the length of the truncated potential range (5σ) plus a small offset,which is 0.2σ in our case. The small offset omits the need to update the cell lists everystep. Here we update the list every 4 steps, which appears to be more often than needed.A cell list stores the current cell for every particle in the simulation, where the particleis located. The neighbor list of cells stores pointers to all neighboring cells of each cell.In fact, it is possible to reduce the neighboring list. This is because we need only paircontributions in the algorithms computing the forces, torques, potential energy, pressureand various correlation functions. Hence, interactions of particles in different cells, whichare already calculated are not considered again in one time step. A further computationspeed increase is achieved by parallelization of the calculations of the pair interactions. Theparallelization and the implementation of cell lists are based on Kohlmeyer (see Appx. D).

C.4. Initial conditions for the molecular dynamics simulations

We start our simulations with particles aligned on a cubic lattice with random Gaussiandistributed velocities and zero angular velocities. In fact, initial velocities and angularvelocities do not matter in a simulation using a thermostat, if the simulations run sufficientlylong to overcome dependency of the initial conditions and this is of course mandatory.The orientations are initially distributed by an uniform random generator. All randomgenerators used are from GSL (see Appx. D).

C.5. Conservation of total momentum in molecular dynamicssimulations

The total momentum in absence of any flow is conserved. However, it is then not possible(or at least not easily possible) to conserve the angular momentum globally, if periodic

115

C. Appendix

boundary conditions are applied. The latter can been understood by assuming a cubicbox with periodic boundary conditions, if a particle leaves the box at one side it appearsat the other side with the same angular momentum. Besides the fact, that the velocityVerlet algorithm (see Sec. 5.1) and the Berendsen thermostat (see Sec. 5.2) in principalconserve the total momentum, it is possible that numerical truncations errors occur andgrow. To ensure numerical accuracy we control the momentum in each spatial directionevery 1000 time step and subtract the Nth part of the numerical error from every velocity,if the absolute value of the total momentum component is bigger than 10−12.

116

DAppendix D.

Programs and libraries

• Libraries: openmp, http://openmp.org/, GSL, http://www.gnu.org/software/gsl/

• Programs: xmgrace, gimp, kile, geany, http://www.gnu.org/

• Programing languages: C++, http://www.gnu.org/

• Script languages: python and vpython, http://python.org/

• Used resources for general ideas of molecular dynamics simulation writing:http://lammps.sandia.gov/index.html (Verlet algorithm, treating dipoles),A. Kohlmeyer optimization-study.pdf, http://sites.google.com/site/akohlmey/software/ljmd(cell list implementation, openmp parallelization)

117

D. Programs and libraries

118

Bibliography

[1] G. Rosenthal and S. H. L. Klapp, J. Chem. Phys. 134, 154707 (2011).

[2] G. Rosenthal, K. E. Gubbins, and S. H. L. Klapp, J. Chem. Phys. 136, 174901 (2012).

[3] I.M. Banat, Bioresource Technology 51, 1 (1995).

[4] G. S. Attard, J. C. Glyde, and C. G. Göltner, Nature (London) 378, 366 (1995).

[5] D. B. Mitzi, Chem. Mater. 13, 3283 (2001).

[6] M. Haumann, H. Koch, P. Hugo, and R. Schomäcker, Appl. Catal. A 225, 239 (2002).

[7] N. Patel and S. A. Egorov, J. Am. Chem. Soc. 127, 14124 (2005).

[8] P. Angelikopoulos and H. Bock, J. Phys. Chem. 112, 13793 (2008).

[9] F. Schmid, Computational Methods in Surface and Colloid Science, edited by M.Borówko (Marcel Dekker, New York, 2000), p. 631.

[10] G. Gompper and M. Schick, Phase Transitions and Critical Phenomena 16, edited byC. Domb and J. L. Lebowitz (Academic, London, 1994).

[11] J. C. Shelley and M. Y. Shelley, Curr. Opin. Colloid Interface Sci. 5, 101 (2000).

[12] H. Bock and K. E. Gubbins, Phys. Rev. Lett. 92, 135701 (2004).

[13] M. Schmidt and C. v. Ferber, Phys. Rev. E 64, 051115 (2001).

[14] J. M. Brader and M. Schmidt, J. Colloid Interface Sci. 281, 495 (2005).

[15] J. M. Brader, C. v. Ferber, and M. Schmidt, Mol. Phys. 101, 2225(2003).

[16] P. S. Christopher and D.W. Oxtoby, J. Chem. Phys. 117, 9502 (2002).

[17] H. Shinto, M. Miyahara, and K. Higashitani, Langmuir 16, 3361 (2000).

[18] R. Goetz and R. Lipowsky, J. Chem. Phys. 108, 7397 (1998).

[19] B. Smit, A. G. Schlijper, L. A.M. Rupert, and N. M. van Os, J. Phys. Chem. 94, 6933(1990).

[20] C. Casagrande, P. Fabre, E. Raphaël, and M. Veyssié, Europhys. Lett. 9, 251 (1989).

[21] P.-G. de Gennes, Angew. Chem., Int. Ed. Engl. 31, 842 (1992).

[22] L. Hong, A. Cacciuto, E. Luijten, and S. Granick, Langmuir 24, 621 (2008).

[23] S. Jiang, Q. Chen,M. Tripathy, E. Luijten, K. S. Schweizer, and S. Granick, Adv.Mater. 22, 1060 (2010).

[24] Q. Chen, J. K. Whitmer, S. Jiang, S. C. Bae, E. Luijten, and S. Granick, Science 331,199 (2011).

[25] F. Sciortino, A. Giacometti, and G. Pastore, Phys. Rev. Lett. 103, 237801 (2009).

119

Bibliography

[26] F. Sciortino, A. Giacometti, and G. Pastore, Phys. Chem. Chem. Phys. 12, 11869(2010).

[27] S. Gangwal, A. Pawar, Ilona K., and O. D. Velev, Soft Matter 6, 1413 (2010).

[28] A. M. Somoza, E. Chacón, L. Mederos, and P. Tarazona, J. Phys.: Condens. Matter7, 5753 (1995).

[29] A. Walther and A. H. E. Müller, Soft Matter 4, 663 (2008).

[30] S. Ye and R. L. Carroll, ACS Appl. Mater. Interfaces 2, 616 (2010).

[31] K. D. Anderson, M. Luo, R. Jakubiak, R. R. Naik, T. J. Bunning, and V. V. Tsukruk,Chem. Mater. 22, 3259 (2010).

[32] A. Perro, S. Reculusa, S. Ravaine, E. Bourgeat-Lami, and E. Duguet, J. Mat. Chem.15, 3745 (2005).

[33] Z. Nie, W. Li, M. Seo, S. Xu, and E. Kumacheva, J. Am. Chem. Soc. 128, (2006).

[34] S. Gangwal, O. J. Cayre1, M Z. Bazant, and O. D. Velev, Phys. Rev. Lett. 100, 058302(2008).

[35] S. Mornet, S. Vasseur, F. Grasset, and E. Duguet, J. Mater. Chem. 14, 2161 (2004).

[36] A. B. Pawar and I. Kretzschmar, Macromol. Rapid Commun. 31, 150 (2010).

[37] D. Suzuki and H. Kawaguchi, Colloid Polym. Sci. 284, 1471 (2006).

[38] S.-H. Kim, S.-J. Jeon, W. C. Jeong, H. S. Park, and S.-M. Yang, Adv. Mater. 20, 4129(2008).

[39] M. D. McConnell, M. J. Kraeutler, S. Yang, and R. J. Composto, Nano Lett. 10, 603(2010).

[40] J. N. Anker and R. Kopelman, Appl. Phys. Lett. 82, 1102 (2003).

[41] C. J. Behrend, J. N. Anker, B. H. McNaughton, M. Brasuel, M. A. Philbert, and R.Kopelman, J. Phys. Chem. B 108, 10408 (2004).

[42] J.-W. Kim, D. Lee, H. C. Shum, and D. A. Weitz, Adv. Mater. 20, 3239 (2008).

[43] N. Glaser, D. J. Adams, A. Böker, and G. Krausch, Langmuir 22, 5227 (2006).

[44] A. Synytska, R. Khanum, L. Ionov, C. Cherif, and C. Bellmann, ACS Appl. Mater.Interfaces 3, 1216 (2011).

[45] J. R. Millman, K. H. Bhatt, B. G. Prevo, and O. D. Velev, Nature Mater. 4, 98 (2005).

[46] V. A. Gilchrist, J. R. Lu, E. Staples, P. Garett, and J. Penfold, Langmuir 15, 250(1999).

[47] R. Steitz, P. Müller-Buschbaum, S. Schemmel, R. Cubitt, and G. H. Findenegg, Euro-phys. Lett. 67, 962 (2004).

[48] F. Tiberg, J. Brinck, and L. M. Grant, Curr. Opin. Colloid Interface Sci. 4, 411 (2000).

[49] O. Dietsch, A. Eltekov, H. Bock, K. E. Gubbins, and G. H. Findenegg, J. Phys. Chem.111, 16045 (2007).

[50] T. Erdmann, M. Kröger, and S. Hess, Phys. Rev. E 67, 041209 (2003).

[51] C. Guerra, A. M. Somoza, and M. M. Telo da Gama, J. Chem. Phys. 109, 1152 (1998).

120

Bibliography

[52] N. Kern and D. Frenkel, J. Chem. Phys. 118, 9882 (2003).

[53] A. Giacometti, F. Lado, J. Largo,G. Pastore, and F. Sciortino, J. Chem. Phys. 131,174114 (2009).

[54] Y. Hirose, S. Komuraa, and Y. Nonomura, J. Chem. Phys. 127, 054707 (2007).

[55] R. Evans, Adv. in Phys 28, 143, (1979).

[56] D. Henderson, Fundamentals of Inhomogeneous Fluids, ISBN 0-8247-8711-0 (MarcelDekker, New York 1992).

[57] Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989).

[58] H. Hansen-Goos and R. Roth, J. Phys.: Condens. Matter 18, 8413 (2006).

[59] R. Roth, J. Phys.: Condens. Matter 22, 063102 (2010).

[60] C. G. Gray and K. E. Gubbins, Theory of Molecular Fluids, ISBN 0-19-855602-0(Clarendon, Oxford 1984).

[61] W. Nolting, Grundkurs Theoretische Physik 6, ISBN 3-540-20505-5 (Springer-Verlag,Berlin Heidelberg, 2005).

[62] R. Steitz, S. Schemmel, H. Shi, and G. H. Findenegg, J. Phys.: Condens. Matter 17,S665 (2005).

[63] P. Hohenberg and W. Kohn, Phys. Rev E 136, B864 (1964).

[64] N. D. Mermin, Phys. Rev. 137, A1441 (1965).

[65] W. F. Saam, C. Ebner, Phys. Rev. A 15, 2566 (1977).

[66] J. A. White, A. González, F. L. Román, and S. Velasco, Phys. Ref. Lett. 84, 1220(2000).

[67] J. A. C. and, Y. Martínez-Ratón, J. Chem. Phys. 107, 6379 (1997).

[68] E. Kierlik and M. L. Rosinberg, Phys. Rev. A 42, 3382 (1990).

[69] J. P. Hansen and I. R. McDonald, Theory of simple liquids, 2nd ed., ISBN 0-12-370535-5 (Academic, London (1986)).

[70] Y. Rosenfeld, M. Schmidt, H. Löwen, 2, and P. Tarazona, Phys. Rev. E 55, 4245(1997).

[71] P. Tarazona, Physica A 306, 243 (2002).

[72] R. Roth, R. Evans, A. Lang, and G. Kahl, J. Phys.: Condens. Matter 14, 12063 (2002).

[73] B. Groh and S. Dietrich, Phys. Rev. Lett. 72, 2422 (1994).

[74] M. Gramzow and S. H. L. Klapp, Phys. Rev E 75, 011605 (2007).

[75] G. M. Range and S. H. L. Klapp, Phys. Rev. E 69, 041201 (2004).

[76] D. C. Rapaport, The Art of Molecular Dynamics Simulation, second edition, ISBN-100-521-82568-7 (Cambridge University Press, 2004).

[77] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, ISBN 0-19-855375-7(Clarendon Press, Oxford 1987).

[78] H.J.C. Berendsen, J. P. M. Postma, W.F. van Gunsteren, A. DiNola, and J.R. Haak,J. Chem. Phys. 81, 3684 (1984).

121

Bibliography

[79] T. Morishita, J. Chem. Phys. 113, 2976 (2000).

[80] D. Frenkel and B. Smit, Understanding Molecular Simulation, ISBN 0-12-267351-4(Academic Press, 2002).

[81] M. Schoen and S. H. L. Klapp, Reviews in computational chemistry 24, ISBN 978-0-470-11281-6 (Wiley-VCH, 2007).

[82] H. Schmidle, C. K. Hall, O. D. Velev, and S. H. L. Klapp, Soft Matter 8, 1521 (2012).

[83] F. W. Starr, J. K. Nielsen, and H. E. Stanley, Phys. Rev. E 62, 579 (2000).

[84] M. Sammalkorpi, S. Sanders, A. Z. Panagiotopoulos, M. Karttunen, and M. Haataja,J. Phys. Chem. B 115, 1403 (2011).

[85] L. Rovigatti, J. Russo, and F. Sciortino, Phys. Rev. Lett. 107, 237801 (2011).

[86] D. S. Rokhsar, N.D. Mermin, and D.C. Wright, Phys. Rev. B. 35, 5487 (1987).

[87] A. L. Mackay, Acta Crystallographica 15, 916 (1962).

[88] J. D. Honeycutt and H. C. Andemen, J. Phys. Chem. 91, 4950 (1987).

[89] J. Farges, M.F. De Feraudy, B. Raoult, and G. Torchet, Surface Science 156, 370(1985).

[90] P. Ilg and E. D. Gado, Soft Matter 7, 163 (2011).

[91] E. Del Gado and W. Kob, Soft Matter 6, 1547 (2010).

[92] R. D. Mountain and D. Thirumalai, J. Chem. Phys. 92, 6116 (1990).

[93] S. H. L. Klapp, Y. Zeng, D. Qu, and R. v. Klitzing, Phys. Rev. Lett. 100, 118303(2008).

[94] S. Klapp and F. Forstmann, Europhys. Lett. 38, 663 (1997).

[95] S. H. L. Klapp and G. N. Patey, J. Chem. Phys. 112, 10949 (2000).

[96] D. Holland-Moritz, D. M. Herlach, and K. Urban, Phys. Rev. Lett. 71, 1196 (1993).

[97] D. Müter, T. Shin, B. Deme, P. Fratzl, O. Paris, and G. H. Findenegg, J. Phys. Chem.Lett. 1, 1442 (2010).

[98] N. R. Tummala and A. Striolo, Phys. Rev. E 80, 021408 (2009).

[99] D. Lugo, J. Oberdisse, M. Karg, R. Schweins, and G. H. Findenegg, Soft Matter 5,2928 (2009).

[100] W. Nolting, Grundkurs Theoretische Physik 4, ISBN 3-540-42116-5 (Springer-Verlag,Berlin Heidelberg, 2003).

[101] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes- The Art of Scientific Computing, 3rd ed. (Cambridge University Press, Cambridge,2007).

[102] M. Ragulskis and L. Ragulskis, Engineering Comput.: Int. J. for Comput-AidedEngineering and Software 23, 368 (2006).

[103] D. Arthur and S. Vassilvitskii, k-means++: The advantages of careful seeding, SODA2007 Proceedings of the eighteenth annual ACM-SIAM symposium, ISBN 978-0-898716-24-5 (2007).

[104] D. L. Davies and D. W. Bouldin, IEEE Trans. Patt. Anal. Machine Intell. PAMI-1,224 (1979).

122

Danksagung

An dieser Stelle möchte ich mich ganz herzlich bei allen bedanken, die mich während meinerDoktorandenzeit mit hilfreichen Diskussionen und natürlich auch moralisch unterstützthaben.

Mein ganz besonderer Dank geht an meine Betreuerin Prof. Dr. Sabine H. L. Klappdie mir meine Promotion ermöglicht und darüber hinaus mich mit vielen konstruktivenGesprächen und Kritiken sehr gut betreut hat.

Ich danke ferner Dr. habil. Thomas Weikl dafür, dass er sich freundlicherweise bereiterklärt hat als Gutachter zu fungieren.

Vielen dank auch an Prof. Dr. Keith E. Gubbins für die freundliche Aufnahme in seineArbeitsgruppe während meines Auslandsaufenthalts in North Carolina (USA).

Für die Organisation und die finanzielle Unterstützung danke ich dem ganzen Interna-tionalen Graduiertenkollege 1524.

Für konstruktive Diskussionen danke ich außerdem Dr. Andrew J. Archer, Prof. Dr. RolandRoth und Prof. Dr. Stefan Sokolowski, sowie während meines Auslandsaufenthalts Lian-gliang Huang, Philipp Kählitz und Jeremy C. Palmer. Hier in Berlin möchte ich ins-besondere Dr. Carlos Alvarez, Thomas Heinemann, Lutz Klaczynski, Nicola Kleppmann,Ken Lichtner, Helge Neitsch und Heiko Schmidle für die hilfreiche fachlichen Diskussionendanken. Ganz besonders danke ich dabei Nicola Kleppmann und Ken Lichtner, die sofreundlich waren und mich mittels Korrekturlesens meiner Arbeit unterstützt haben. DesWeiteren danke ich auch Simon Dinter der sich ebenfalls ein paar Kapitel durchgelesen hat.Darüber hinaus geht mein Dank an unsere ganze Arbeitsgruppe für ein sehr freundlichesUmfeld und nette drei Jahre. Zu guter Letzt möchte ich mich noch bei meiner Familiebedanken, die mir immer motivierend beiseite stand.

Theory and computer simulations of amphiphilic Janus particles · 2017-10-26 · We apply molecular...

Documents

Transcript of Theory and computer simulations of amphiphilic Janus particles · 2017-10-26 · We apply molecular...