Molecular Dynamics simulations of the morphology of

1 downloads 0 Views 1MB Size Report
3.3 The Flory Huggins mean field approximation . . . . . . . . . . . 22 ... avenue to tailor-make materials. Experiment plays a central role in ...... Society A: Mathematical, Physical and Engineering Sciences, 368, p. 1009. 8. Filomia, F.; Faginas Lago, ...
Molecular Dynamics simulations of the morphology of mixtures of organic molecules and/or polymers.

Chemistry Colloquium by Rodrigo Perez-Garcia Supervised by

Dr. Alex H. de Vries (RuG)

Abstract

The aim of this report is to provide an overview of the role of Molecular Dynamic (MD) simulations as applied to mixtures of polymers and carbon based compounds. The text is divided in four parts.

The rst part is a contextualization of

the present situation regarding morphology simulations for the kinds of mixtures proposed. The second part focuses on pointing out some of the dierent modeling techniques available. The third part of the report aims to explain the position that MD occupies within this framework and the relationship with the other techniques. A set of case studies will be investigated more in depth in the fourth part. The selection of this set of case studies was selected with a focus on aspects of personal relevance for future work, recent publication-time and expected applicability. The last part of the report provides a series of conclusions on the topic of MD advantages and disadvantages.

Contents

1 2

3

4

Introduction

4

Modeling techniques

8

2.1

Quantum techniques . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2

Classical molecular and atomistic methods . . . . . . . . . . . . .

11

2.2.1

Historical overview . . . . . . . . . . . . . . . . . . . . . .

11

2.2.2

Molecular Dynamics . . . . . . . . . . . . . . . . . . . . .

12

2.2.3

Types of ensembles . . . . . . . . . . . . . . . . . . . . . .

14

2.2.4

Molecular mechanics force eld . . . . . . . . . . . . . . .

14

2.2.4.1

All-atom MD . . . . . . . . . . . . . . . . . . . .

14

2.2.4.2

Coarse grain MD . . . . . . . . . . . . . . . . . .

15

2.3

Microscale methods . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.4

Mesoscale and macroscale methods. . . . . . . . . . . . . . . . . .

18

Molecular Dynamics as a key middle-point: Modeling mixtures 20 3.1

Polymer and organic molecules mixtures . . . . . . . . . . . . . .

22

3.2

Modeling a system

. . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.3

The FloryHuggins mean eld approximation . . . . . . . . . . .

22

Case studies 4.1

25

A molecular simulations study of the miscibility in binary mixtures of polymers and low molecular weight molecules. S. S. Patnaik, R. Pachter, (Polymer, 2002)

4.2

. . . . . . . . . . . . . . . . .

25

Characterization of polymer-fullerene mixtures for organic photovoltaics by systematically coarse-grained molecular simulations. D. M. Huang, A. J. Moule, R. Faller, (Fluid Phase Equilibria, 2011).

4.3

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Molecular dynamics simulation study of binary fullerene mixtures. R. Ruberto, M. C. Abramo, C. Caccamo, (Physical Review, 2004)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

5

Conclusions

36

6

References

38

2

Foreword

To PGF, GPF, JPG, PPG, BGM, CFR and CBR

According to Heidegger, philosophy arises from (and returns to) the existence.

Philosophy is thus the very analysis of existence.

This work speaks

of a part of this factual existence of mine, thrown into the world of life so to speak. Science and philosophy lead the path of knowledge in which more or less successfully I am trying to walk through. Truth, as stated in various dictionaries say, is a property of some of our ideas. This means the latter adapt to reality. Falseness means unsuitability of our ideas to reality.

Scientic research as I understand it, in an broad and nontrivial

sense, can be described as the search for demonstrable truths. This is somehow in contradiction with the scientic ideal, with the myth, sometimes almost a religious belief, that science is about discovering nal, infallible, absolute or ultimate truth. That myth, if openly endorsed by scientists, inadvertently or otherwise, is dangerous because it fuels gross conicts like impact-level-factor inconsistencies and increases public confusion about science. It is important nowadays to t theoretical predictions to observables in order to reduce costs and to enhance possible research discoveries, coherently drawing a line (as well as this can be done) between statements, or systems of statements, of the empirical sciences and the theoretical disciplines. What I mean is to create models that can, as accurate as possible, describe Reality in silico as a this will be the way to achieve solutions faster and more eciently. This, together with open access to publicly nanced knowledge and synergistic collaborations, will lead us to a brighter future in as broad a meaning as possible. What I mean is to improve eciency and solving active problems of current research remembering that The whole is more than the sum of its parts (Aristotle).

3

Chapter 1 Introduction

A major breakthrough in chemistry was the postulation and subsequent discovery of macromolecules by the German chemist Hermann Staudinger in the rst quarter of the 19th century.

In 1953, German chemist Karl Ziegler

patented the method for polymerization of ethylene (gure 1.1), and in 1954, the Italian Giulio Natta developed polypropylene (gure 1.2), which are the two most widely used plastics today. In 1963, these two scientists shared the Nobel Prize in Chemistry for their studies on polymers.

The role of polymers and

organic molecules has gained importance with time and are nowadays widely used in everyday life [1,2].

Polymer research constitutes a broad research eld, consisting of topics such as: the selection of a solvent for use with a solute, control of solutes dissolution, diusion and swelling processes, the particular modiers and their use and the study of phase diagrams of molecular and polymeric mixtures. Similar to polymers research, the design of other organic molecules for certain precise uses, it provides a very powerful tool in many dierent

elds of chemistry. The devel-

opment of composite polymer-small organic molecule systems provides a further avenue to tailor-make materials. Experiment plays a central role in science. It is the wealth of experimental results that provides a basis for understanding the chemical behavior of these

Polyethylene fragment. Ziegler was able to polymerize ethylene under ambient conditions using trialkyl aluminum and titanium chloride as catalyst. Figure 1.1:

4

Polypropylene fragment. Natta developed polypropylene, a isotactic polymer. It is a strong, lightweight, and heat-resistant plastic. The material is durable but can be molded into thin objects that bends without breaking.

Figure 1.2:

systems we are focusing on. Experimental techniques, such as X-ray diraction or nuclear magnetic resonance (NMR), allow determination of the structures and elucidation of the interplay of each component of the mixtures. Yet, experiment is possible only in conjunction with theories. Computer simulations have altered the interplay between experiment and theory.

The essence of the simulation

is the use of computers to model a physical system.

Testing properties of a

system using computer modeling is faster and less expensive than synthesizing and characterizing it by laboratory experiments. It is mandatory to bear in mind there are dierent problems arising as molecules exhibit a wide range of time scales over which specic processes occur; for example: - Local Motions (0.01 to 5Å,

10−15 to 10−1 s)

which include Atomic uctua-

tions, Sidechain Motions, Loop Motions. - Rigid Body Motions (1 to 10Å,

10−9 to 1s) Helix Motions, Domain Motions

(hinge bending), Subunit motions. - Large-Scale Motions (> 5Å,

10−7 to 104

s) Helix coil transitions, Dissocia-

tion/Association, Folding and Unfolding. There are also properties such as electronic conductivity, photochemistry, chemical degradation that that can not be described or explained even qualitatively without considering the electronic wavefunction.

On the other hand,

other properties (polarizability, molecular charge distribution, etc.) are usually treated eectively averaging over the electronic degrees of freedom. Bulk properties, are dependent on various factors mainly due to the arrangement of the molecules and the packing and this needs longer simulation times. Each one of these processes should be tackled with dierent levels of theory. When treating mixtures bulk properties need to be studied and the calculation becomes more expensive.

Hence, the applicability of the techniques available

appear to be restricted (gure 1.3). If the goal of reliably predicting properties from rst principles via atomistic simulations, using intermolecular potentials [3] derived from quantum chemistry, could be turned into reality, many industries would prot, and clearly this would also be a major success of statistical mechanics. A thorough theoretical assessment of the properties of these binary systems is important, but dicult due to the lack of availability of accurate parameterization for MD simulations [4].

It is also important to bear in mind that

obtaining appropriate results for the required time- and length scales of the systems is dicult using most available methods.

The problem is even more

complicated when one of the components is polymeric because additional factors such as molecular packing, chain- exibility and molecular weight need to be considered [5]. Since 1982 experimental and theoretical advances have been made by sev-

Figure 1.3:

ology

Graph presenting the correlation of structure under stody vs. necessary method-

eral groups [6,7,8,9] in both the area of performing simulations of mixtures and through making various accurate force elds available for Molecular Dynamics (MD). It is also important to note that micro- and mesoscale methods are bringing some useful insight into the behavior of binary carbon-based systems' microstructure. However, the theoretical understanding of compressible binary uid mixtures is a dicult problem even for mixtures of medium-size molecules [10] and up to now it has not been done using QM/MD methodologies. Quantum Mechanical (QM) electronic structure calculations on mixtures of relatively small molecules have been performed by several researchers [11,12], while for complex fullerene mixtures it has yet to be examined.

However, in

ref.[13], the authors claim to have successfully performed dynamics of the interaction of fullerene and one hydrogen molecule (gure 1.4). Despite the successful implementation of this approach, a major problem is that by using ab

initio methods of a very high degree of sophistication, the calculation becomes extremely time consuming and the estimate of the thermodynamic properties of other correlated mixtures cannot yet be accomplished on a routine basis. The principal problem is therefore, as already stated, the diculty of establishing adequate forms of the intermolecular potentials [6], since already a single polymer coil exhibits nontrivial structure from the scale of a chemical bond to the coil radius (100 Å) which is already complex to model, and the simulation of collective phenomena required huge simulation boxes containing many polymers, and as a result calculations are therefore very expensive. Although much work has been done studying mixtures, nothing can be found that deals perfectly with these systems. Therefore, this literature review will try to enumerate what has been done, giving a general overview that will help in the placement of the problem of mixture processes in an appropriate context. The basic elements of the dierent techniques useful for describing mixtures of polymers and organic molecules will be introduced, to, then, arrive to the basic principles of molecular simulation techniques (Chapter 2), mainly molecular dynamics (Chapter 3) - going through the possible interactions, force elds available, ensembles and Flory mixing thermodynamics. Quantum techniques

The 3D H2 probability distribution of (a) one and (b) two p-H2 inside C70, from quantum diusion Monte Carlo calculations [13] Figure 1.4:

and microscale methodologies are also presented (Chapter 2), albeit less detailed. Finally, a series of case studies (Chapter 4) are reviewed, that will help to understand the potential uses of the available techniques in this eld, in the hope it will serve as attention enhancer and perspective broadener rather than as source of confusion.

To nally give a series conclusions (Chapter 5) and

references to the literature revised (Chapter 6).

Chapter 2 Modeling techniques

Polymeric and organic structures in melts, blends or mixtures can have sizes that range from nanometers to microns, millimeters and larger. The corresponding time scales of the dynamic processes relevant to unravel dierent properties span an even wider range, from femtoseconds to milliseconds or even seconds or hours for glassy materials or for large scale ordering processes such as phase separation in blends.

No model or simulation algorithm can span this range

of length and time scales.

Therefore, molecular and mesoscopic models for

polymeric materials range from those including quantum eects and electronic degrees of freedom; to chemically realistic, classical models; to coarse-grained, particle-based mesoscale models that retain only the most essential elements of the polymer system to be simulated; to

eld-theoretic models that describe

the polymer system in terms of density or composition variables (see

gure

2.1). One solution, which holds particular challenges for polymer materials, is multiscale simulation, the bridging of length and time scales and the linking of computational methods to predict macroscopic properties and behavior from fundamental molecular processes. Thus, based on accurate QM calculations, a force eld is determined, which includes charges, force constants, polarization, van der Waals interactions and other quantities that accurately reproduce the QM calculations.

With the FF, the molecular dynamics (MD) are described

with Newton's equations of motion instead of the time dependent Schrödinger Equation.

The MD level allows predicting the structures and properties for

much larger systems (in terms of number of atoms) than for QM, allowing direct simulations for the properties of many interesting systems. This leads to many relevant and useful results in mixture behavior; however, many critical problems in this eld still require time and length scales far too large for practical MD. Hence, there is a need to model the system at the mesoscale (a scale between the atomistic and the macroscopic) and to pass messages from the atomistic scale to the mesoscale and to even the macroscale. This linking through the mesoscale in which the microstructure can be described is probably the greatest challenge to developing reliable rst principles methods for practical blends description. Only by establishing this connection from microscale to mesoscale it is possible to build rst principles methods for describing the properties of mixtures and composites. Any solution for passing information from one scale to another (upper) scale

8

Figure 2.1: Molecular modeling and simulation methods commonly used for polymer nanocomposites, adapted from [14].

is based on the denition of multiscale modeling which only consider 'objects' that are relevant at that particular scale, implicitly describing selected degrees of freedom of smaller scales by summarizing those degrees of freedom by using some representative parameters. All MD approaches are initially based on the application of a force eld that transfers information on the dierent energy terms that describe the deviation of bond lengths, bond angles and torsion angles away from equilibrium values, plus terms for non-bonded pairs of atoms describing van der Waals and electrostatic interactions from quantum chemistry to atomistic simulation. From atomistic simulation to a macroscale one (a traditional approach based on the estimation of for example the characteristic ratio, the Kuhn length, and the Flory Huggins interaction parameter, can be used) to explore the macroscopic properties of a system through these microscopic simulations, for example, to calculate changes in the binding free energy of a particular polymer candidate, or to examine the energetics and mechanisms of conformational change. The connection between microscopic simulations and macroscopic properties is made via statistical mechanics which provides the rigorous mathematical expressions that relate macroscopic properties to the distribution and motion of the atoms and molecules. This approach for determining the input parameters for large mesoscale simulations is based on the following information: (i) the bead size and chain architecture, (ii) the bead mobility M, and (iii) the eective Flory-Huggins

(χ)

interaction parameter.

2.1

Quantum techniques

Computational techniques that provide more condence and rigor are those based on quantum mechanics (QM). With the development of computational power and dierent techniques of numerical calculus, computational modeling has shown its capabilities as an eective way study a great variety of molecules and their behavior with the so called ab initio molecular dynamics method. Several ways can be used to predict dynamic behavior of chemical systems

at an electronic level.

QM is able to generate an accurate description of the

potential energy surfaces on which to propagate then the dynamics. Quantum computational methods were developed to solve the electronic Schrödinger equation that results from the time-independent Schrödinger equation (see eq.2.1), normally working under frozen nuclei approximation (or more formally Born-Oppenheimer approximation).

ˆ = Eψ Hψ

(2.1)

In practice there are three main methodologies used to obtain the potential energy surface (E of eq. 2.1 as a function of the coordinates) where to propagate the dynamic: ab initio, semiempirical and those based on Density Functional Theory (DFT). The philosophy behind each of the three methods is quite dierent. Ab initio methods use quantum chemical formalism rigorously, using no further assumptions in addition to the rst principles (the postulates of quantum mechanics). The method of self-consistent eld or Hartree-Fock (HF), is the simplest among the ab initio. It tries to nd the wave function that minimizes the energy of the system. Based on the variational theorem (the energy obtained is

E0

ε ≥ E0 , where

is the ground state of the Hamiltonian), it is possible to calculate the wave

function considering that each electron moves in a eld due to the nuclei and there is an average eld created by all other electrons as a monodeterminant orbital function. Thus, the HF method does not explicitly include the interaction between each electron and the rest of the other single electrons (electronic correlation) but an average potential, which causes a overestimation of the electron repulsion.

The simplest way to include the electronic correlation is the

conguration interaction (CI). This represents the electronic ground state of the molecule as a mixture of all possible electronic states (Full-CI calculations). However, in practice the expansion must be truncated to a nite (small) number of congurations (truncated-CI calculations) otherwise the time required for these calculations would be extremely long. In general, CI is often considered a benchmark method, and often provides one of the best non-relativistic descriptions possible for a molecular system. expensive for medium size systems.

However, this method is prohibitively

Another way to treat electronic correla-

tion methods is perturbation theory. Undoubtedly, the most widely used is the Møller- Plesset (MP), which is based on the theory of Rayleigh and Schrödinger. These procedures are, however, computationally time-consuming, and can not deal with big-systems dynamics [10,11,12,15]. Semiempirical methods were born with the purpose of providing quantummechanical techniques for the study of molecular properties, being suciently accurate and reliable to have a practical value in chemical research, without the disadvantages of ab initio procedures, i.e. they apply to larger molecules or molecular systems. These methods introduce a number of simplications in the expression of the HF Hamiltonian, by the introduction of a set of parameters (like for example assuming certain Dierential Diatomic Overlaps). Among the most popular semiempirical methods available at the present time it is possible to highlight the Austin Model 1 (AM1) [16] and the Parameterization Method 3 (PM3) [17] developed by Dewar and Stewart, respectively. The DFT approach is based on the theorem of Hohenberg-Kohn. It states that the monoelectronic density can be used to obtain all the observables (for instance E) of the ground state of a system. This is empirically obtained using the

Kohn-Sham method [18]. In addition, the electron density can be represented as the sum of

N

one-electron orbital densities, which makes DFT calculations

easily applicable [14]. In 1985, Quantum Molecular Dynamics (QMD) simulations appear, pioneered by Car and Parrinello, under the DFT framework [19]. This technique, as stated by Shubin Liu, still is probably the most eective technique to compute the dynamic properties of materials from electronic structure equations with a reasonable computer time [20]. It calculates the nucleus classically, take a classical step in agreement with Newton's laws and then calculates the electronic structure that provides data for the new calculation of the nuclear positions, again classically. It is a technique that calculate the dynamics on the y. It is important to notice that QM methods are able to describe both ground states and excited states, where, in the latter, no experiment could be used to fully understand the behavior. However, the methods with less rigor, such as semiempirical, in my opinion, will slowly disappear in favor of next level simulation methods.

2.2

Classical molecular and atomistic methods

2.2.1 Historical overview The simulation by molecular dynamics is a technique used to calculate the equilibrium properties and transport of a classical many-body system. It was rst introduced by Alder and Wainwright in the late 1950's [21]. Given an N-body system which has no analytical solution, the only way of solving it is by using a numerical solution method. The algorithm of a molecular dynamics program consists of the numerical solution of the equations of motion of the particles under consideration. The solution of the equations provides a realistic path (coordinates and conjugate momenta as a function of time) of the system under study. An important issue is the relationship between the properties of matter whether in solid, liquid or gas - and the interaction between atoms or molecules that constitute this material. Instead of trying to deduce the microscopic behavior and predict directly from the experiment, the molecular dynamics method tries to reproduce the behavior using parametrized models. Classical Molecular Dynamics methods are based on the following principles:

ˆ ˆ

Nuclei and electrons are treated collectively as atom-like particles These atoms are spherical-like particles (sizes are obtained from measure-

ments or theory) and can have a partial charge appropriate for the liquid state (possibly obtained from theory).

ˆ

The bonded particles interact as "springs" and the interactions are repre-

sented by classical potentials.

ˆ ˆ

These interactions must be predened for specic sets of atoms The interactions determine the spatial distribution of particles and their

energy. These models consider the molecular chemical dynamics of atoms as spherical masses and bonds as properly parametrized springs. The mathematics of spring deformation can be used to describe the ability of the links to "stretch"

The atoms in a molecule can be modeled as charged spheres connected by springs that maintain bond lengths and angles. The charged atoms interact with each other (via Coulomb's law) and also with the particular solvent. The model shown is called a molecular mechanics potential energy function. The molecule shown is an alanine dipeptide. From www.stanford.edu/~boas/science/protein_design/index.html. Figure 2.2:

"bond twist" and "dihedral twist". Non-bonded (more than two bonds apart) atoms, dovetail through van der Waals forces, steric repulsion and electrostatic interactions. These non-bonding interactions are easiest to describe mathematically when atoms are considered as spheres with characteristic radii. Molecular chemical dynamics will therefore predict the energy associated with a given conformation of a molecule based in these interactions.

However, the energy

of molecular dynamics are not absolute quantities. Only energy dierences between two or more conformations have a physical meaning. The equation for the total energy of a given system under study together with the data (parameters) needed to describe the behavior of dierent types of atoms and bonds, is called "force eld". Here (gure 2.2) it is shown that a standard force eld contains bonded terms relating to atoms that are linked by covalent bonds, and nonbonded terms describing the long-range electrostatic and Van der Waals forces.

2.2.2 Molecular Dynamics Classical MD models do not take into account all the details of the electronic structure at quantum level and thus cannot predict the inuence of specic groups of atoms on the physical properties as accurately. However, in the last few years, atomic scale MD has been quite successful in exploring dierent phenomena occurring at pico- to nanosecond time scales [22,23]. The structure and behavior of mixtures at the atomistic level can be simulated via two distinct methods, the pure Molecular Dynamics (MD) and the Monte Carlo (MC) simulation methods.

In an MD simulation, the motion of individual atoms

within an assembly of N atoms or molecules is modeled on the basis of either a Newtonian deterministic dynamics or for example a Langevin-type stochastic dynamics, given the initial position coordinates and velocities of the atoms [24]. In a common MD simulation, a nite model of a mixture consisting of

N

atoms

and conned in a simulation cell of volume

V

is initially constructed. The cell

is replicated in all spatial dimensions generating its own periodic images as well as those of the original

N

atoms.

The imposition of this periodic boundary

condition is necessary in order to compensate for the undesirable eects of the articial surfaces associated with the nite size of the simulated system. The energetics and dynamics of the atoms are obtained from inferred interatomic potentials,

φ(rij ),

and the simulation involves the computation of forces expe-

rienced by each atom at each simulation time step (see table 2.1) from these potentials. A common guideline when simulating exible molecules is to choose a time step at least one order of magnitude smaller than the time of the shortest period of motion.

System under study

Types of motion

Time step

Atoms Rigid Molecules Flexible Molecules (rigid bonds) Flexible Molecules (exible bonds)

Translation Translation, Rotation Translation, Rotation, Torsion Translation, Rotation, Torsion, Vibration

10−14

Table 2.1:

5x10−15 2x10−15 −15

10

to5x10−16

Suggested time steps (of the dierent types of motion) for various systems [14].

In a MD simulation involving a complex mixture whose atoms interact via a complicated interatomic potentials, the simplifying assumption (it will also help to save computation time) that each atom interacts with its nearest neighbors, located in both its own cell and the image cells that lie within a specied cut-o radius, is usually made. The interacting

N

atoms require

3N

simultane-

ous coupled second-order dierential equations of motion. These equations are integrated numerically within the computational cell by a variety of numerical integration techniques, all based on the nite dierence methods, where time is discretized. MD, at this level, allows one to investigate into the eects of mixing dierent materials and/or the dynamics close to the interface [14]. The Monte Carlo technique is a stochastic method in which the structural facts and the thermodynamic properties are obtained from averages.

Thus,

the Monte Carlo method attempts to obtain a Boltzmann sampling of a system according to an algorithm that consists of generating random structures (link with microscale Brownian Dynamics method), deciding on the basis of an appropriate energy criterion whether the new structure is accepted or disregarded. The simpler algorithm formalism is called Monte Carlo Metropolis, in which the generation of new structures takes place at random, ensuring microscopic reversibility of the system [14]. However, in dense chemical systems, such as polymers, the Metropolis method is inecient and one must use more sophisticated algorithms. In this sense the algorithm called Congurational Bias, which ignores the existing barriers and performs a random walk over local energy minima is a good approach to calculate geometries and properties of the phase molecular distributions [14,25]. In general, both Molecular dynamics and Monte Carlo methods have become important tools for the study of uids and mixtures. They have been used to study the equilibrium and transport properties of model atomic blends and models for molecular liquids.

2.2.3 Types of ensembles An ensemble is a group of all possible systems which have dierent microscopic states but have an identical macroscopic or nal thermodynamic state. There are dierent ensembles with dierent characteristics.

ˆ Microcanonical ensemble (NVE) : The thermodynamic state characterized N , a xed volume, V , and a xed energy, E . This

by a xed number of atoms,

corresponds to an isolated system.

ˆ Canonical Ensemble (NVT): This is a collection of all systems whose therN , a xed volume, V , and a xed temperature, T .

modynamic state is characterized by a xed number of atoms,

ˆ

Isobaric-Isothermal Ensemble (NPT): This ensemble is characterized by a

xed number of atoms,

ˆ

N,

a xed pressure,

P,

and a xed temperature,

T.

Grand canonical Ensemble (µVT): The thermodynamic state for this en-

semble is characterized by a xed chemical potential, a xed temperature,

µ,

a xed volume,

V,

and

T.

2.2.4 Molecular mechanics force eld A molecular mechanics force eld [25,26,27] analytically describes the potential energy of a system in terms of the geometries of atomic centers.

The

energy calculation and dynamics force eld can be derived from ab initio quantum mechanics and empirically (using spectroscopy and crystallography). We will explain in more detail all-atom and coarse grain force elds (see gure 2.3).

2.2.4.1

All-atom MD

In all-atom MD, all of the atoms (including hydrogens) are treated explicitly during the calculations. It is the most realistic approach, and generally prevails over united atoms (e.g., methyl groups treated as a single units [28]), implicit solvent (e.g., distance dependent dielectric or other approximations to account for the lack of solvent molecules), and other methods with reduced complexity. More recently, increases in computer speed, and the proliferation of inexpensive multi-processor machines have enabled all-atom simulations of full proteins on long (a few microseconds) simulation time scales [28,29]. All atom simulation techniques provide atomic resolution of equilibrium protein dynamic behavior and non-equilibrium events. When used in conjunction with experimental disciplines, simulations provide an enhanced view of the system under study.

Visual comparison of an all-atom structure (left) and a coarse grain description of the same (right), of alamethicin peptide. In order to describe the coarse grain system, MARTINI CG force-eld was used. The force-eld has been adequately parameterized [39,40] in agreement to a truthful representation of the chemical structure. The backbone parameters (back-bone bonded terms) are dependent on the secondary structure of the beads but independent of the individual amino acid. Secondary structure was forced by including a dihedral potential between backbone atoms.

Figure 2.3:

2.2.4.2

Coarse grain MD

Together with atomistic simulations, for certain processes, coarse-grained [27-28] models allow exploration of larger size systems and longer timescales at the same computational cost.

As the fast and slow molecular dynamics

are suciently independent, it should be possible to ignore fast vibrations for the study of slow dynamics.

Thus, to study the mechanisms, dynamics and

structural changes of several systems, as those systems presented here, some level of coarse grained (CG) description of the system can be applied as it is shown in the second case study. CG allows to study dynamical phenomena in a mixture by adjusting parameters such that the bonded potentials, diusivity and self-diusion coecient data, of the same system all-atom MD, are reproduced as close as possible [28].

The representation of all protein amino acids and other molecules mapped into beads (Figure from http://md.chem.rug.nl/cgmartini/). Figure 2.4:

More precisely, in CG models, molecules are described by the interaction sites representing groups of atoms, providing a reduced resolution description of a given system.

These models are expected to be highly computationally

ecient, both because mapping atoms into the sites reduce some of the degrees of freedom, and also because high frequency intramolecular vibrations have been incorporated into averaged eective interactions between sites. Consequently, it is possible to choose a larger time step∆t, and therefore speed up computations as stated above. There are many dierent levels of CG models, ranging from the 'united atoms' approach (CH2 and

CH3

groups are considered as a single

unit), to mesoscale models using the rigid regions of well-dened equilibrium structures previously analyzed with higher level of theory. This broad range is one of the assets of MD. An important level of system description is the  residue resolution level represented by the MARTINI CG model [27], where the atoms of each residue are mapped into beads (gure 2.4). It was initially developed for lipids in water simulations and later optimized and validated for proteins [28] and also very recently for mixtures of nanocompounds and dierent solvents [29]. MARTINI CG has become very popular due to its accurate parametrization of a large library of biologically relevant building blocks. In CG models with  residue level resolution, the lack of proper hydrogen bond description can make it necessary to combine CG model with diverse methods to maintain the overall shape of the system under study. Although the dynamics of the atomistic model is deterministic and conservative, the dynamics of the resulting coarse-grained model is stochastic and includes dissipation and possibly memory eects [30]. However, using Coarse grain models is becoming of paramount importance as it substantially decreases the number of the system's degrees of freedom, provoking an increase in the velocity of the simulations by several orders of magnitude [28].

2.3

Microscale methods

The next (downgrade) level of theory will be a structure evolution analysis which implies microsecond simulation times.

Microscale methods are the

connection point between molecular and continuum regimes. Some mechanisms and properties must be interrogated using atomistic and related concurrent multiscale methods, many others may be best examined using mesoscale or microscale methods. They incorporate the molecular description implicitly and describe the kinetics with a respectable level of exactness. Within this frame there are various methodologies like Brownian Dynamics, Dissipative Particle Dynamics (DPD), Latice Boltzman (LB), Time dependent Ginzburg-Landau Method (TDGL) or Dynamic DFT method. Brownian dynamics theory is born within the idea of Brownian motion (g. 2.5) and random generation of possible structures. It neglects inertia and the state of the polymer is therefore completely specied by the positions of the

N +1

beads,

ri .

Hydrodynamic interactions between the beads are introduced

in a pairwise-additive approximation through the mobility matrix connects the mean velocity (eq. 2.2) of bead

i

to the force on bead

µij , j,

which

Figure 2.5: Example of Brownian motion generated using Mathematica software. x =0 ; y = 0; brownian2D = Table[{x = x + ( 2 *Random[] - 1 ) * d r , y = y + ( 2 *Random[] - 1 ) * d r ) , { i f 1000)]

vi =

X

µij · Fj

(2.2)

Fj = −∇rj φ The conservative force,

Fj

(2.3)

(see eq. 2.3), is derived from the potential energy

of the polymer, that consists of

N +1

beads connected by springs between

neighboring beads. The dynamics of a polymer in a low-molecular weight melt can be described using this ingeniously simple model, where the beads move as Brownian particles subject to random forces and to frictional forces proportional to their velocity exerted from their environment. DPD, which has exactly the same structure as conventional Brownian motion dynamics, with the important dierence between the two dynamics that the Brownian dynamics does not satisfy Newton's third law and therefore does not preserve time track-back and conservation of overall momentum.

The latter

is the main drawback of these microscale methods (not reproducibility), since the dynamic Brownian particle velocity are coupled directly to the dissipative force[14].

DPD describes diusion properties well but fails for hydrodynamic

ow. Lattice Boltzman is mainly used in polymers solution dynamics.

The LB

methods were developed to study the dynamics of complex mixtures. In contrast to particle-based methods such as Dissipative Particle Dynamics, the uid information is stored on a regular grid, simplifying the computation in a parallel programming environment. In common with other methods that include the uid degrees of freedom explicitly, LB simulations are very exible and can include Brownian motion via thermally-driven uctuations applied directly to the uid.

Polymer solutions are modeled by coupling the polymer to an LB

uid with a frictional drag force [7]; this method has been applied to conned polymers as well. The time dependent Ginzburg-Landau (TDGL) method is a microscale method used to simulate the structural evolution of phase separation.

Scale comparative. At molecular lever, atomic conformation is appreciated. At the microscale, domains are observed. At the mesoscale, grains with numerous domains. At the macroscale of the specimen, is possible to observe numerous grains with dierent domain structures. Figure 2.6:

2.4

Mesoscale and macroscale methods.

Mesoscale modeling uses a basic unit just above the molecular scale, and is particularly useful for studying the behavior of polymers and soft materials. It can model even larger molecular systems, but with the commensurate trade-o in accuracy (gure 2.6).

Furthermore, it is possible to transfer the

simulated mesoscopic structure to nite elements modeling tools for calculating macroscopic properties for the systems of interest. These methods obey the basic laws of continuity (conservation of mass), equilibrium, momentum principle and conservation of energy and entropy [14]. Continuum mechanics deals with ideal homogeneous mixtures.

Its aim is

to describe their response to external inputs using appropriate constitutive relations. Both mesoscale and macroscale methods aim at nding a volume element's responses to prescribed mechanical forces and to deduce from them the corresponding overall properties. The most straightforward application of such studies is materials characterization. In order to get realistic phase distributions, the analysis of the spatial variations of the microelds (in suciently large reference volumes) tends to be beyond present capabilities, and approximations have to be used.

For conve-

nience, the majority of the resulting modeling approaches may be treated as falling into two groups: The rst of them comprises methods that describe the microgeometries of inhomogeneous materials on the basis of statistical information:

ˆ

Mean Field Approaches (MFAs) and related methods: Such descriptions

typically use information on the microscopic topology, inclusion shape and orientation, and also on the statistics of the phase distribution. Mean eld approaches tend to be formulated in terms of the phase concentration tensors, they pose relatively low computational requirements, they have been highly successful in describing the elastic response of inhomogeneous materials.

ˆ

Variational Bounding Methods: Variational principles are used to obtain

upper and lower bounds on the overall elastic tensors, elastic moduli, secant moduli, and other physical properties of inhomogeneous mixtures. Many analytical bounds are obtained on the basis of phase-wise constant stress (polarization) elds. Many bounding methods are closely related to MFAs. Because none of them explicitly account for n-particle interactions, Mean eld approaches are sometimes referred to as  noninteracting approximations in the literature. They postulate the existence of a proper representative vol-

ume element (RVE) and typically assume some idealized statistics of the phase arrangement at the microscale. The second group of approximations is based on studying discrete microgeometries aiming at fully accounting for the interactions between phases.

It

includes [14]:

ˆ

Periodic Microeld Approaches (PMAs), also referred to as periodic ho-

mogenization schemes or unit cell methods.

In these methods the inhomoge-

neous material is approximated by an innitely extended model material with a periodic phase arrangement.

The resulting periodic microelds are usually

evaluated by analyzing repeating unit cells providing a sort of detailed overview of the mixture, but they tend to be computationally expensive.

ˆ

Embedded Cell or Embedding Approaches (ECAs): A model consisting of

a  core containing a discrete phase arrangement that is embedded within some outer region to which force eld loads are applied. The material properties of this outer region may be described by some macroscopic constitutive law, they can be determined self-consistently from the behavior of the core, or the embedding region may take the form of a coarse description and/or discretization of the phase arrangement.

ECAs can be used for materials characterization,

and they are usually the best choice for studying regions of special interest in inhomogeneous materials.

ˆ

Windowing Approaches: Typically rectangular or hexahedral shape sub-

regions ( windows ) are randomly chosen from a given phase arrangement and subjected to boundary conditions that guarantee energy equivalence between the micro- and macroscales. Further descriptions of inhomogeneous systems are not discussed here due to their typically rather weak physical background and their rather limited predictive capabilities of mixtures behavior.

Chapter 3 Molecular Dynamics as a key middle-point: Modeling mixtures

Multiscale modeling can be dened as the enabling technology of science and engineering that links phenomena, models, and information between various scales of complex systems [14,22,23,24]. Molecular Dynamics (both all-atom and Coarse Grain) play an important role in this multiscale linking. The idea of using MD modeling in multiscaling is straightforward: one computes information at a smaller (ner) scale and passes it to a model at a larger (coarser) scale omitting degrees of freedom. Thus, based on accurate Force Field (FF) parametrization, which includes charges, force constants, polarization, van der Waals interactions and other quantities that accurately reproduce the QM calculations, MD can reproduce/perform accurate calculations in less time. Using a FF, the dynamics of the studied melt or mixture is described according to Newton's equations, instead of the Schrödinger Equation. The MD level allows predicting the structures and properties for systems much larger in terms of number of atoms than for QM, allowing direct simulations for the properties of many interesting systems.

This leads to many relevant and useful results

in materials design; however, many critical problems in this eld still require time and length scales far too large for practical MD. Hence, there is a need to model the system at the mesoscale (a scale between the atomistic and the macroscopic) that links the atomistic scale to the mesoscale and the macroscale. This linking through the mesoscale in which the microstructure can be described is probably the greatest challenge to developing reliable rst principles methods for practical application of materials' design. Only by establishing this connection from microscale to mesoscale it is possible to build rst principles methods for describing the properties of new materials.

The problem here is that the

methods of coarsening the description from atomistic to mesoscale or mesoscale to continuum are not as obvious as it is in going from electrons to atoms. In other words, the coarsening from QM to MD relies on basic principles and can be easily generalized in a method and in a procedure, while the coarsening at higher scales is system specic [23].

20

These methods employ a potential energy function (referred to as a force eld) that calculates the overall potential energy of the system based on the summation of individual atom-atom pair interactions. The force eld equation takes into account the contributions due to bonded interactions (bond stretching, bending, and torsion) as well as non-bonded interactions (van der Waals and electrostatic). These energy contributions are determined by a set of empirical parameters that are used by the force eld to calculate the energy contribution for each type of interaction for the atom types that are dened in the molecular system under consideration. To accurately simulate molecular behavior using an empirical force eld, parameters of the force eld must be balanced and tuned together to appropriately represent the behavior of the given molecular system being addressed. As a consequence of this, a force eld designed for one type of application cannot be condently applied to other applications without separate validation; this issue is known as force eld transferability. QM to MD approaches are initially based on the application of a Force Field that is able to reproduce properties from quantum chemistry in the atomistic simulation. From atomistic simulation to mesoscale one can use a traditional approach based on the estimation of the characteristic ratio of the mean square end-to-end distance (the Kuhn length), and the Flory Huggins interaction parameter.

This approach for determining the input parameters for mesoscale

simulation is based on the following information: (a) the bead size and chain architecture, (b) the bead mobility M, and (c) the eective Flory-Huggins X interaction parameters. The molecular dynamics equations of motion, for all but the simplest systems, are of sucient complexity that the integration must be done numerically over a large number of very small discrete time steps rather than analytically. This treatment of time assumes that at any given discrete time step the atomic forces are xed. This assumption holds only if the magnitude of the time step is suciently small. At any given time step, these coordinates are used to calculate the potential energy and its rst derivative, the force, using a molecular mechanics force eld. Generally, for any atom, evolution in time proceeds from step n to n + 1 as described in equations 3.1 to 3.3, where subscripts denote the time step,

f

∆t is the magnitude of the integration time step, a is the acceleration, m is the atomic mass, v is the velocity, and x refers

is the force on the atom,

to the atomic coordinates:

an =

f(x,y,z) m

(3.1)

vn+1 = vn + an ∆t

(3.2)

1 xn+1 = xn + vn ∆t + an ∆t2 2

(3.3)

A long series of these steps generates a trajectory through phase space, the

6N

dimensional space (where

N

is the number of atoms) dened by the three

space vectors of the atom positions and velocities. In general, post-simulation analysis regards the atomic position (coordinates) subspace of the phase space.

3.1

Polymer and organic molecules mixtures

The dynamics of a large polymer molecule are complex.

This complexity

also happens for organic molecules of big size. A dynamics study of a mixture of both components is, therefore, non trivial.

Even though individual atoms

move approximately with the same equilibrium distribution of speeds as if they were somehow independent, their motion is constrained by the existing chemical bonds. As the structure of the components gets bigger it needs longer times to rearrange in a stable phase distribution, and quantities depending on overall conformation, such as overall distances or bulk conformation, take a very long time to change from their original values. Because of the long relaxation times associated with large-scale conformational rearrangements, enhanced by the fact that polymeric uids have memory associated to the liquid state macroscopic recoil. Therefore, physically meaningful results, comparable to experiment, about structure (diusion coecient, viscosity, etc,) together with the thermodynamic properties, given the interaction potentials between the atoms, imply using models capable of analyzing specic interactions that occur at time scales comparable with the relaxation times of the whole chain.

It is important to quantitative

understand properties, processingperformance and structure relations in materials. On the other hand, in order to reproduce motion correctly one has to solve the equations of motion. The timestep has to be chosen as a compromise between fast simulation and satisfying the equilibrium condition.

3.2

Modeling a system

The main ingredient of a simulation is a physical model. For a molecular dynamics simulation, this means choosing a potential or force eld, function

V (r1 , ...rN )

of the positions of the nuclei, representing the potential energy of

the system when the atoms are arranged in a specic conguration. The forces are then derived as the gradient of the potential with respect to atomic movements as it is shown in eq. 3.4:

Fj = -∇ri

XX i

φ(ri − rj )

The simplest choice for V is written as a sum of pair interactions The second summation index

j > i,

(3.4)

j>i

φ(ri − rj ).

tells us to consider each pair just once.

The choice of force eld depends in large part on the system to be studied and its properties to be investigated. Beyond rened description of the relevant atomic vibrations, the conformations of macromolecules depend more on the twists, the repulsion, van der Waals attractions and electrostatic interactions.

3.3

The FloryHuggins mean eld approximation

Thermodynamics of polymeric systems with dissimilarities between dierent molecules in a blend or mixture, can be well studied using the Flory-Huggins

theory.

The Flory-Huggins theory of polymer solutions is based on a lattice

theory. In 1941 Flory and Huggins, independently, developed a simple theory able to explain non-ideal behavior of polymer mixtures. Flory made a complex derivation but got a very simple and intuitive result for polymer-solvent and polymer-polymer systems [31], eqs. 3.6 and 3.7. Given the lattice volume fractions

φ1

φ1 = N1 /N

φ2 ,

eq. 3.5

φ2 = xN2 /N

(3.5)

∆Hm = RT n1 φ2 χ12

(3.6)

∆Gm = RT [ n1 ln φ1 + n2 ln φ2 · x1 + n1 φ2 χ12 ]

(3.7)

N1 is the number of N2 is the number of polymer molecules, φ2 is the volume fraction of the polymer, x1 is the number of segments in the solvent molecule, n indicates the number of moles and R is the gas constant. The value of the interaction parameter comes from the solubility parameters, δi , and the molar volume of the solvent, v , eq. 3.8. Where

χ12

and

is the Flory-Huggins interaction parameter,

solvent molecules and

2

χ12 =

v (δ1 − δ2 ) RT

(3.8)

Their model is based on a thought experiment since it is based on the occupations of a lattice where molecules can not be given an exact position. The Flory-Huggins interaction parameter, eq. 3.8, lets us understand the energetics in the mixing of polymers because it indicates the strength of unfavorable interactions. A phase diagram can be constructed to summarize the phase behavior of the mixture, showing regions of stability, instability, and metastability calculated from Flory-Huggins. This theory is useful when trying to investigate the eects of

χ12

on molecular weight, and content of each component on the

microscopic morphologies, mixing kinetics, and viscosity of blends. Therefore, it is very useful for MD as results show that Flory-Huggins interaction parameters computed by atomistic MD simulations are consistent with experiment and little computationally demanding.

Composition ratio, temperature, and

Flory-Huggins are optimal parameters particularly the latter bridge very well the gap from atomistic simulations, where solubility parameters can be calculated, to mesoscopic simulations where mesophases and network formation can be studied. In the mean eld approximation the interactions between molecules are assumed to be due to the interaction of a given molecule and an average eld produced by all the other molecules in the system.

To aid in modeling, the

solution is imagined to be divided into a set of cells within which molecules or parts of molecules can be placed (lattice theory, see gure 2.1).

In order to generate a low enthalpy of mixing, the dierence in solubility parameters must be close to zero and the molar volume must be very small. The smaller the molar volume, the larger the solvent molecule. However, the larger the molecule (larger values for mixing expansion.

n1

and

x1 ),

the bigger the enthalpy of

In order to maximize solubility, a polymer solvent must

either be very small, or have an equivalent solubility parameter to the polymer. Although this lattice model was originally developed for polymers, it can be extended to a include mixtures in a more general way. As stated above, this theory provides approximate relations for the heat of mixing and the entropy. From these energies several thermo-statistical quantities can be thus derived, like activities and number density uctuations. Since, in short, the interaction parameter measures the anity of a solvent towards a determined polymer it is a very relevant parameter when studying mixing. There are several examples where Flory-Huggins is applied successfully when trying to explain a mixture's properties such heat and volumes of mixing, lower critical solution temperatures, and the enthalpic and entropic components of the chemical potential. One of the most cited examples are by Lacombe and Sanchez [32, 33] where the authors show how benecial it would be to use the FH theory (and its modication) to predict the previously mentioned properties. As this theory is simpler than other methodologies and it is therefore widely used to get this values for the coarser MD simulation.

Chapter 4 Case studies

The analysis of these case studies intends to bring an understanding of (increasingly) complex systems. We will try to determine and dene the research questions presented in each one of the papers. Then, we will add strength to what is already known through previous research by evaluating the applicability and the advantages/disadvantages of each approach. In each of the the case studies, the target is to emphasize a detailed contextual analysis of a limited number of conditions (depending on the paper) and their relationships. To some extent we will try to portray a complex problem in a way that conveys the state

of the art of it.

4.1

A molecular simulations study of the miscibility in binary mixtures of polymers and low molecular weight molecules.

Pachter, (Polymer, 2002)

S. S. Patnaik, R.

Flory-Huggins parameter is seldom used for describing miscibility. The free energy of mixing

∆Gmix

and the interaction parameter

χ

are very important

in determining whether a binary system is miscible. A negative mixing energy indicates that the mixing is exothermic and that the two components tend to mix together. In contrast, a positive mixing energy means that the two components tend to separate. Similarly, if If

χ

χ

is negative, the two components prefer to mix.

is positive, the two components prefer to phase separate.

In this study it is presented how the Flory-Huggins parameter can be obtained from Molecular Dynamics simulations.

It is important to be able to

model it due to the implications the parameter has towards the understanding of mixtures and these system's many dierent applications. It is good to remark this parameter is of paramount importance to gain insight into the processes of miscibility of small molecules and polymers. The article presents an analysis of a mixture of poly(methylmethacrylate): PMMA, and the liquid crystal (4-n-pentyl-4'-cyano biphenyl): 5CB (see gure

25

(a) Skeletal structure of methyl methacrylate, the monomer of PMMA

constituent

Figure 4.1:

4.1).

(b) Skeletal structure of 4-n-pentil.4'-cyano biphenyl

Schematic representation of the molecules used

This is a blend of great importance for the understanding of polymer

dispersed liquid crystals, PDLC. It has been found that the morphology of the polymerization process is highly dependent on where in the phase region it is initiated [34].

In order

to study this process, and due to the importance of an accurate description of these binary mixtures, dierent simulation methods are used for the analysis, mainly the thermodynamic phase behavior by Flory-Huggins and Molecular simulations.

To investigate these and other related questions with computer

simulations require models rich enough to capture:

the intrinsic free energy,

the main structural features, and the connection within the structures and the interactions between them. The model must also be simple enough to allow for an ecient simulation of several molecules in the solution, which limits the use of atomistically detailed descriptions. The negative free energy (see eq. 4.1&4.2) of mixing can be used to predict the mixing as a function of composition for dierent temperatures.

∆Hmix = ∆Emix + P ∆Vmix

(4.1)

∆Gmix = ∆Hmix − T ∆S

(4.2)

It is also presented in the article how joining two theories, the previously mention Flory-Huggins (FH) free energy method [31] and the Maier-Saupe (MS) free energy of nematic ordering (microscopic model for phase transitions applicable for liquid crystals) [35], is necessary to describe the miscibility process. In the case of isotropic mixing, MS can be neglected and once the FH parameter is known, the entire phase behavior of the system can be calculated. Flory and Huggins introduced the idea of dividing the free energy, or chemical potential of mixing, into a combinatorial term and with a residual (contact or interchange) term which depends only on the local environment of each unit of polymer. The combinatorial term is associated with the number of ways of placing the molecules into the mixed environment. Therefore the FH can describe mixing assuming the polymer is formed by a series of connected segments occupying a regular lattice.

It is important to notice that, here, the interac-

tion parameter initially proposed by Flory and Huggins is replaced with a more general one, depending on the temperature, pressure and concentration.

Due to the dierences in size of the system constituents, it is complicated to accurately dene the topologies in the presented system [37]. In the paper it is stated that many advances have taken place by providing accurate force elds (correct parametrization), which is useful as it provides insight into the behavior of binary systems [37]. It is important to acknowledge this study as it does not want to be an absolutely complete basis, but it is a valid path to a better understanding of atomic interactions. Therefore although MD can not be used routinely to predict phase behavior, it can certainly provide a better fundamental understanding of miscibility. They have obtained the interaction parameter (ignoring concentration dependency) of 5CB with PMMA and methyl methacrylate (both monomer and oligomer) using two dierent simulation techniques in order to evaluate the validity / applicability of each one.

The interaction parameter in Eq.

3.8 can

be derived using dierent simulation methods, most of which are based on the assumption that the molecules interact through pair interaction energies. The rst method, called the two segment method, considers the molecules as isolated pairs of molecules and calculates, using a MC approach, a broad range of orientations (coordination numbers, torsional angles,...) and hence an average interaction energy. The calculation is carried out using Blends module in

CERIUS2

and the COMPASS force eld. The dierential energy, calculated

using FH models for each pair, adds up to obtain

∆Gmix

following the Eq.

4.3. This is a fast way to evaluate miscibility without neglecting many specic interactions.

∆Em = RT φ1 φ2 χ12

(4.3)

As stated earlier the Monte Carlo - MD method is appropriate because random walk over local energy minima is a good approach to calculate average geometries and properties for mixtures. However, as stated in Patnaik's paper, bulk eects are poorly described with MC. Therefore, it is not recommended to determine the solubility using this method (mainly for the bigger systems presented) as

χ

(depending in the coordination number) would be meaningful

only when the conformation of the shorter polymeric segments are representative of the big polymer. Cohesive energy densities (CED) were the parameters evaluated in the second place. CED in liquid mixtures usually depend on the chemical structures, especially on those constituent groups of the molecules present. CED are dened in a way they show the energy necessary to break all intramolecular bonds in a volume-unit and hence quantify the relative strengths of the interactions. Then these cohesive energies derived from MD allow the determination of the interaction parameter, assuming no specic interactions are formed nor destroyed. Co1

hesive energy (Ecoh ) can be expressed (using Hildebrand relationship in terms of a solubility parameter

δ.

The

Ecoh

2 ) δ = E coh

for pure compounds or mixtures

is calculated from the congurations obtained during the MD simulation after equilibration was reached. The

Ecoh

can be calculated by subtracting the to-

tal conguration energy of the pairs in one phase from the total conguration energy of one mole of pairs in the other phase, which is directly related with the eq. 3.8. In this way the conguration energy corresponding to the internal binding and conformational energies of the molecules (from bonds, angles,

Figure 4.2:

Calculated cohesive energy vs. molecular weight

dihedral angles, and internal non-bonded interactions) will not be taken into account.

There is plenty of experimental data for

δ

and this makes it more

simple to verify the computed interaction parameter, which makes this method increasingly applicable. In this second case the program used is

DISCOVERs always

with COM-

PASS force eld under NTV conditions (further description of the simulation conditions is provided in their paper). They agree that comparing individual cohesive energies to estimate miscibility is not as good as calculating the energy of mixing between dierent systems from variations of the cohesive energy. The second approach in their paper further improved their initial methodology. Calculating the cohesive energies, deriving them at dierent temperatures and concentrations for both mixed and unmixed systems. It is much more computationally demanding but shows the dependence of

χ

with respect to the two

parameters. The simulation set up is similar to the previous approach, and the results obtained concur with experimental values. However, the improvements in the analysis of the miscibility are not high, although some extra insight is provided for discovered concentration dependence of the interaction parameter. Next picture (gure 4.2) is showing how miscibility decrease as the monomer polymerize. The cohesive energies can be derived at dierent temperatures and concentrations and both temperature and concentration dependence of

χ

can

be studied knowing that:

∆Emix = −(Ecoh )AB + φA (Ecoh )A + φB (Ecoh )B Summarizing, rstly, using the two segment method to analyze local interaction energies achieves a fairly good insight on the dierent contributions to the mixture. The diculty in dening a single coordination number for binary mixtures limits the applicability of this method. Consequently, it was developed an idea to use modied coordination numbers that depend both on the relative size of the segments and the concentration of the mixture. Although this is a computationally fast method, the nal results were found to be very dependent on the choice of polymer fragment, concluding that a short polymeric segment was not representative of the PMMA conformation. The second approach manages to calculate cohesive energies, with and without taking into account the dependence of

χ

with the concentration and tem-

perature. The demanding calculations were found not to be extremely advantageous in comparison with those of the previous method. Their conclusion is that with the advances in atomistic molecular dynamics simulations and with better force elds, it is possible to describe miscibility. Even if the atomistic simulations provide the most reliable results, convening in an adequate way the available simulating techniques will provide a good insight in progressively-increasing-size polymer miscibility.

This parametrization can

be used in mesoscale models to study morphology evolution in time.

Conclusions Results reported in this work show preliminary conclusions obtained in the rst stages of a wide computational eld like for instance mixtures of polymers and small organic molecules. By carrying out these molecular simulation approaches for the prediction of thermo-physical properties of multicomponent mixtures, Patnaik and Patcher give a new detailed description of the process that can be used in the future. The rst studies in this eld are those of ref. [37]. They help the reader by placing in context molecular simulation methods for calculation of phase equilibrium from assumed intermolecular interactions. There has been rapid progress in the eld over the past decade. Before the mid-1980's, even the calculation of a phase diagram for a potential as simple as the Lennard-Jones potential was considered a very complex task. The current paper shows a study, at dierent levels of theory, on the importance of dierent eects that must be taken into account to describe atomic level interactions of miscibility suitable for future calculations. This is meaningful only when the conformation of the shorter polymeric segment is representative of the polymer. We know from comparison with experimental results that properties extrapolation from the conformation of the short polymeric segment (in isolation) is not representative of PMMA. This was indeed conrmed by further analysis. The torsional angles were extracted from MD simulations with periodic boundary conditions such as the analysis for determining the miscibility of methyl methacrylate monomer/5CB was reasonable successful using two-segment method, but not so for PMMA/5CB. Subsequently a dierent method that allows for bulk eects to be included was considered. This is indeed one of the strengths of MD. The validity of the methodology is baked-up with the many experimental data from the cohesive energy of the dierent organic molecules and polymers used. This makes this study quite valuable for the scientic community and reliable for predicting bulk properties of polymeric systems. Up to now it is been cited 16 times in diverse topics ranging from verication of FH applicability [39], copolymerization [39,40], miscibility of blends [40,41] and even for determine anti-carcinogenic drug's solubility [42].

4.2

Characterization of polymer-fullerene mixtures for organic photovoltaics by systematically coarse-

D. M. Huang, A. J. Moule, R. Faller, (Fluid Phase Equilibria, 2011). grained molecular simulations.

In this work a computational framework is developed to study qualitative morphology evolution of mixtures close to those used in Organic Photovoltaics (OPV) (gure 4.3). OPV morphology and their theoretical inference is a key element that has been studied thoroughly [25,43]. Characterization of the morphology of these devices is a challenging scientic problem that few techniques can address. Experimental evaporation rates, solvent type, etc, are not easy to visualize and quantify, therefore these long simulations (mostly achievable by CG-MD) are a powerful tool which could provide the required spatial resolution sensitive to morphological details present at bulk as well as the composition. The feasibility of Coarse Graining this kind of mixtures and its evolution towards bulk-heterojunction formation is outlined here.

Snapshots of the simulation conguration of the system with P3HT:C60 = 1.27:1 (w/w) with Nmono = 48 at t = 0 ns (left), t = 30 ns (center), and t = 135 ns (right). The C60 molecules in the largest cluster are highlighted in blue and all other particles in the system are shown as dots. The length of each side of the simulation box is roughly 25 nm. Taken from [38]. Figure 4.3:

CG-MD allow various process and variables to be controlled independently, as well giving an automated exploration of the phase space. It will be helpful for tailoring manufacturing process as it gives insight in many of the applicable properties and at the same time nal arrangements can be used for future quantum calculations. Study of phase separation in a system of P3HT/C60 of several exciton diusion lengths in size (i.e. is approaching the device scale) require a not small computational eort that CG-MD can handle with easiness. Modeling this system using adequate coarse grain methods will help to move towards a more complete understanding of phase and morphology particularities of these increasingly important systems. In order to prepare an organic photovoltaic device (OPV) a disordered bulk-heterojunction needs, usually, to be generated.

Most common polymer-

based solar cells are based in the mixture of the conducting polymer (poly(3hexylthiophene): P3HT) and substituted fullerene structures. In this particular paper by Huang et al., a systematic coarse grain (CG) MD model for the rather big conjugated polymer - fullerene systems is developed. They, however, do not

add functionality to the fullerene as the chemical intuition may point out as the idea is not exatly model a useful system but to set the basics for a sistematic approach to these systems. Experimentally, these structures are obtained by drying a mixture of the two species, previously casted, onto a substrate.

The arising network of do-

mains, separate (~10 nm) creating the appropriate structure able to generate a photocurrent after light absorption [44]. OPV are easy to produce from P3HT and PCBM (fullerene derivate) mixtures in a volatile solvent.

Dierent mor-

phologies are possible. The dierent morphologies strongly aect to the nal eciency of OPV devices. The inability to experimentally visualize morphology evolution hinders the ability to quantify the eect of various process and system variables (such as evaporation rate, blend ratio, solvent type) on simulated morphology evolution. Therefore, there is a big interest to understand and predict the structure of these systems from rst principles. There are computational approaches to this problem, but are mostly limited to one scale: atomistic scale or macroscale. From a macroscale perspective the problem of thin lm formation is well-studied, as it is summarized by Norman et al.

[26].

They explain the macroscale lm thickening during evaporation

process with angular velocity of the coater, concentration, evaporation rate, and solution viscosity. However, morphology evolution is not analyzed in this study. This, as proposed here, can be studied using MD as a computational framework that can model morphology evolution at the nanoscale and at the same time represent macroscale phenomena. For experimentally relevant system's sizes, an all-atom simulation would be very computationally expensive, hence the use of MD in general and CG systems in particular became conditio sine qua non while describing a system with 100,000 sites (corresponding to ~1.2 million atoms).

Figure 4.4: Chemical structures of P3HT and C60 with coarse-grain transformation (in red). Taken from [38]. The study focuses on morphology prediction not evaluating electronic structure as it intends to be a systematic approach to these mixtures. This particular system does not have much applicability, as conrmed by prof. Hummelen in a personal communication, due to the high tendency of C60 to crystallize. By using the Iterative Boltzmann Inversion (IBI) methodology the whole

(a) ∝Radius of gyration (solid lines) and ∝End-

to-End distances (dashed lines) as a function of time, starting from a pseudo-random conguration of polymer chains and C60 molecules, for various P3HT:C60 blend ratios and polymer chain lengths: P3HT:C60 = 1:0 (w/w) with Nmono = 24 (crosses) and Nmono = 48 (circles), P3HT:C60 = 1.85:1 (w/w) with N = 48 (squares), and P3HT:C60 = 1.27:1 (w/w) with Nmono = 12 (up triangles), Nmono = 24 (down triangles), and Nmono = 48 (diamonds). The temperature T is 550 K in all cases. Figure 4.5:

(b) Probability distribution p(κ2) why of polymer chain shape anisotropy for various P3HT:C60 blend ratios and polymer chain lengths: P3HT:C60 = 1:0 (w/w) with Nmono = 48 (circles) and P3HT:C60 = 1.27:1 (w/w) with Nmono = 12 (up triangles), Nmono = 24 (down triangles), and Nmono = 48 (diamonds). The temperature is 550 K. (Only distributions for systems where equilibration could be ensured are shown.) Taken from [34].

system was systematically coarse grained (see gure 4.4) from atomistic simulations by the same group [43]. All calculations were performed using a massively parallel calculus software (Lammps MD package), running under NPT conditions. The IBI approach start with a test simulation. It is run with an arbitrary chosen initial potential.

This results in a trial radial distribution functions

(RDF) for this rst optimization step. Taking the ratio between trial and target RDF leads via Boltzmann-inversion to a potential correction.

By adding

the correction to the old potential, the new potential for the next iteration is generated.

Eective intra- and intermolecular potentials are then rened

through an iterative process, which eventually causes the distribution and correlation functions from coarse-grained simulation to converge to their atomistic ally computed counterparts. It is important to remark their previous study [43] shows a great agreement of their CG and the atomistic description. The parameters under study are the size of the chains in the melt (see g 4.5(a)) and the shape anisotropy of chains (see g 4.5(b)) (calculated using the average radius of gyration and end-to-end distances).

From (a) it is possible

to conclude that the polymers chains are separated from fullerenes. The latter weakly inuence the overall polymeric structure. The information from (b) is that long chains assure anisotropy, hence equilibration.

This can be inferred

knowing that broad range of conformation are visited for high values of Nmono. Although the simulations performed are too short to observe full phase separation, local fullerene cluster formation is perceived and studied. The clusters have a big surface area. This is good as continuous paths within both phases are needed to be able to establish a photocurrent. It is important to notice the time scales for atomistic and CG simulations dier (ve to one respectively). Thus, although the CG time scales in our simulations are not directly comparable

with the  real time scales, previous results [43] suggest that these time scales are likely of the same order of magnitude. Even if the present study can only describe phase separation qualitatively, a further development in the quality of the model should lead to a quantitative description.

Conclusions Their model can describe the local morphology in bulk heterojunctions on the length scales of several tens of nanometers, close to the photovoltaically relevant scale.

The validity of this coarse grained method is no doubt [43].

These large-scale simulations will in the future be able to describe, following their coarse grained method, the morphology over substituted fullerenes. This work is the continuation of ref. [43] and in conclusion, the simplied model described in this work constitutes one of the rsts attempts to satisfactory deal with the dynamics of polymer fullerene mixtures and is protable both for the authors of this paper that continue to work on similar mixtures [45] and for other four working groups dealing with bulk-heterojunction modeling. I believe interesting to point out that in some of these four the study of Huang et al. is mention as an model that bring insight in bulk heterojunction simulation [46,47].

4.3

Molecular dynamics simulation study of bi-

R. Ruberto, M. C. Abramo, C. Caccamo, (Physical Review, 2004)

nary fullerene mixtures.

Here it is presented a preliminary study about the phase behavior of supercooled and solid solutions of dierent size fullerenes in order to analyze in quantitative detail the structure and dynamic evolution of the mixture. Taking into account that the experimental studies regarding the substitutional solid solution of C60 and C70 are controversial, various calculations in

silico have been done in order to achieve more denitive results. Particularly a study of C60-C70 mixture (at 50%) using MC based on LJ potential representation [48], predict solid solution. The Girifalco model [49] that is currently used for fullerenes and carbon nanotubes [56], describes the interaction in terms of analytically spherically symmetric pair potentials. This allows to map C60 and its phase diagrams in a study published in 2003 by Cheng et al. Cases with

Cj

(j=70, 76, 84, 96) have been examined more recently. Within this work, adopting the conditions of Girifalco, it is intended to perform extensive constant pressure MD of C60-C70 mixtures for evaluating on them the solid phase and miscibility gaps, among other properties. The study has later been extended to other three dierent size ratios diameters

(α)

depending in the

σi . α=

σC60 ; (j = 70, 76, 84, 96) σCj

If we think on C60 as good approximation of an isotropic sphere (symmetry

Ih )

with diameter 7.1Å, the simplest, and perhaps more intuitive, way to treat

these systems is using the well studied Hard Spheres (HS) model.

The other

Radial distribution functions (rdf) at room temperature of C60-C70 solid solutions, at dierent concentrations. (Rdfs are sifted by two units on x-axis) [25]. Figure 4.6:

Cj

are easily linked with spherical shapes also.

The steeply repulsive nature

of Girifalco potential make it possible to treat these systems following the HS model. Mapping of the MD calculated phase behavior of fullerene mixtures are found to evolve similarly to HS systems presenting some insolubility regions, in agreements with ref. [50]. It is important to mention that, in this paper, pairwise molecular interactions were described by the Girifalco potential in ref. [49] presented in a slightly dierent form (Kniaz et al. [51]) to accommodate attractions and repulsions between spheres of dierent size and charge. It reads as follows:

υ

=



+

Ni

Ni Nj A 48di dj r

"

Ni Nj B 36di dj r

"

1 [r − (di + dj )]3 1 [r − (di + dj )]9





1 [r + (di − dj )]3 1 [r + (di − dj )]9





1 [r − (di − dj )]3 1 [r − (di − dj )]9

stand for the number of C atoms and and

di

+

+

#

1 [r + (di + dj )]3

#

1 [r + (di + dj )]9

ith υ=

for the radius of the

A and B are the parameters entering the 12-6 potential B . The parameters A and B are those determined by Girifalco et al. for 12 r r pure C60. A and B will remain the same for all Cj , the diameter di is assumed fullerene species.

=6A +

to scale with the number of carbon atoms. As mentioned before four dierent

α,

at eleven concentration proportions are analyzed using NPT conditions. In their study, crystallization is analyzed. It takes place by forming a highly

defective fcc structure with lattice sites occupied at random by the two fullerene species. They evaluate parameters such as crystallization temperature, RT density, volume, Enthalpy, to show how these evolve in correlation to the concentration. Two main results are reached, regions where specic values for the concentration are particularly favored (C60 and C70 rich ends) and a range between

x ≤ 0.5 where a very dierent behavior is shown,

0.3 ≤

probably related with the fact

of amorphous phases being formed (see gure 4.6). The patterns for C60-C70 mixtures have a similar shape to an eutectic coexistence line as in a HS model. A correspondence of the eutectic point, in fairly good agreement with HS mode, shows the existence of colloidal glasses. For the other values ofα the region of insolubility extends through a broader range tending to zero for the C60-C96 mixture. This lower boundary is predicted also by Home-Rothery rule [52]. These emerging similarities with colloidal-like behavior could be explained, as stated in a paper of 2004 by the same authors, due to sort and/or macroparticle interactions in the complex system.

Conclusions The eect of the size dierences between fullerene species is big: The solid in-miscibility region rapidly expands even upon a tiny reduction of

α,

with

formation of an amorphous phase at relatively low temperatures (~1100K). Experimental studies of these systems still encounter some diculties (very dicult to control

α

ratio, extreme conditions, ...)

and therefore computer

simulation methods have become valuable tools for exploratory research on these systems. Despite the usefulness of the understanding of the behavior of fullerene mixtures, little theoretical work had been done before the publication of this paper. The organic materials deriving from fullerenes are of high importance in new materials design.

These mixtures are an almost ideal toy-system for

investigations of phase behavior, and testing the predictions of other theoretical models and computational experiments. It is therefore important to notice that fullerenes appear to be more soluble when the two components of the mixture have similar sizes in any range of concentrations except for a middle-range where the fullerenes arrange themselves in a close packed liquid phase.

This study provide a realistic picture of the

morphology and phase structure that is consistent with both observations and other models and may help explain the behavior of this increasingly important systems. The simplicity and the computational economy of the HS model will allow for these simulations a much larger system over a longer period of time at a low computational cost.

Chapter 5 Conclusions

The advantage of MD approach is that the detailed level of description allows deeper study of systems' behavior unachievable using other methods. Even though modeling has already had great success and holds even greater promises in a wide variety of areas, many challenges remain in order to make it into a solid scientic discipline perfectly coupled with experimental data. The usefulness of QM varies considerably depending on what sort of properties are to be calculated.

The same occurs for MD and semi-macroscopic

methods. It is therefore very intelligent to use multiscale modeling to be able to cover the gaps from atomistic to the macroscopic scales in reasonable times. Mathematical study of the multiscale models and consistency with the full microscale model is of great importance.

To assure adaptivity by developing

new style of error indicators to guide the renement algorithms and to check which parameters are easily transferable from one level to another. As we could see in one of the case studies, Flory-Huggins parameter works very well and is easily used in transferring information of mixtures. The IBI coarse graining methodology also allows to systematically connect atomistic and CG methods and is a good referential point. The existence of microscales models (as HS) must be taken into account as they can simplify the manifold and dicult challenges arising in the simulation of certain mixtures properties.

These rather drastic

approximation does lead to paradigmatic results. The main limitation of molecular dynamics is, at present time, the lack of intermolecular potentials that can describe adequately, for example, complex mixtures. A priori quantum mechanical techniques are not yet in position to give accurate enough potentials for calculations in the kind of systems proposed due to size restrictions so using all-atom or coarse-grain MD is of paramount importance as these parameters can be obtained from other sources and the results are in good agreement with reality. Even though the most promising approach appears to be using multi-scale modeling with parameters tted to experimental phase equilibrium and thermodynamic data available, further validation must be done. It will allow the linking between theoretical levels. Microscale methods are starting to be applied when studying continuous systems. These are promising as they allow accurate determinations of critical points which can be, therefore, used for higher level technique. Therefore, combining the results of molecular simulations as middle point, continuum modeling

36

and quantum simulations will increment the performance and the accuracy. The nal aim is to develop computational techniques that both accurately predict materials properties and provide an estimate of their accuracy, and then to apply these methods to accelerated materials development and enhance the understanding of the eect of atomic-scale processes on meso and macroscale behavior.

MD fulll this perfectly and the development progression is rather

high which open a window of hope for these technique. As a summary of the dierent levels of theory and their applicability see table 5.1:

Quantum Dynamics (relativistic and non-) ⇐⇒Atomic Nuclei, Electrons, Photons Quantum Dynamics (atomic) ⇐⇒Atoms, Ions, Molecules (photons) LEVEL c: Molecular Dynamics⇐⇒ (Macro-)Molecules, Fluids, Solutions, Liquid Crystals LEVEL d: Brownian Dynamics⇐⇒ Non-Equilibrium Processes, Colloidal systems, Polymers LEVEL e: Fluid Dynamics⇐⇒Non-Equilibrium Macroscopic Fluids, Gases, Solids

LEVEL a:

LEVEL b:

This table shows how reduction of complexity decreases as the level of theory increases, adapted from H.J.C. Berendsen [53].

Table 5.1:

Chapter 6 References

1. Paul J. Flory - Nobel Lecture December 11, 1974. Nobelprize.org 2. Karl Ziegler - Nobel Lecture December 12, 1963. Nobelprize.org 3. Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham W. A.; (1981) Intermolecular Forces, Their Origin and Determination. Clarendon Press, Oxford. 4. Patnaik, S. S.; Pachter, R.; (2002) A molecular simulations study of

the miscibility in binary mixtures of polymers and low molecular weight molecules Polymer, 43, p. 415. 5. Olabisi, O.; Robeson, L.M.; Shaw, M.T. (1979) Polymer-polymer miscibility New York: Academic Press. 6. Swinton F. L.; (1987) Mixtures of non-polar molecules Pure &AppI. Chem., Vol. 59, No. 1, p. 35 7. Higgins, J. S.; Lipson, J. E. G.; White, R. P.; (2010) A simple approach to polymer mixture miscibility Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 368, p. 1009. 8. Filomia, F.; Faginas Lago, N.; A Simplied Myoglobin Model for Molecular Dynamics Calculations Computational Science and Its Applications - ICCSA 2006, Springer Berlin / Heidelberg, 3980, p. 731. 9. Perlmutter, J. D.; Drasler II, W. J.; Xie, W.; Gao, J.; Popot, J-L.; Sachs, J, N.; (2011) All-Atom and Coarse-Grained Molecular Dynamics Simulations of a Membrane Protein Stabilizing Polymer Langmuir, 27 (17), p. 10523 10. Rowlinson, J. S.; Swinton, F. L.; Liquids and Liquid Mixtures, Butterworths, London 1982. 11. Marcon, V.; Breiby, D. W.; Pisula, W.; Dahl, J.; Kirkpatrick, J.; Patwardhan, S.; Grozema, F.; Andrienko, D. (2009) Understanding structuremobility relations for perylene tetracarboxdiimide derivatives. Chem. Soc., 131, p. 11426

38

J. Am.

12. Cheung, D. L.; McMahon, D. P.; Troisi, A. (2009). Computational study of structure and charge transfer parameters of low-molecular mass P3HT. J. Phys. Chem. B, 113, p. 9393 13. Sebastianelli, F.; Xu, M.; Ba£i¢, Z.; Lawler, R.;Turro, N.J.; 'Hydrogen Molecules inside Fullerene C70: Quantum Dynamics, Energetics, Maximum Occupancy, And Comparison with C60' (2010) Journal of the American Chemical Society, 132 (28), p. 9826 14. Zeng, Q.H.; Yu, A.B.; Lu, G.Q., (2008) Multiscale modeling and simulation of polymer nanocomposites.

Progress in Polymer Science, 33, p.

2 15. Irle, S.; Zheng, G.; Elstner, M.; Morokuma, K.; (2003) Formation of Fullerene Molecules from Carbon Nanotubes: A Quantum Chemical Molecular Dynamics Study. Nano Lett., 3, p. 465 16. Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P.; 'AM1 : A New General Purpose Quantum Mechanical Molecular Model' (1985) Journal of the American Chemical Society, 107, p. 3902 17. Stewart, J. J. P.; (1989) Optimization of parameters for semiempirical methods I. Method J. Comput. Chem. 10, p. 209 18. Kohn, W.; Sham, L. J.; (1965) Self-Consistent Equations Including Exchange and Correlation Eects Physical Review, 140 (4A), p. 1133 19. Car, R.; Parrinello, M.; (1985) Unied Approach for Molecular Dynamics and Density-Functional Theory Phys. Rev. Lett. 55, p. 2471 20. www.sccas.cn/gb/cooper/cooper16/reports/Shubin.ppt 21. web.mit.edu/mbuehler/www/.../lecture_1_nano-macro-buehler.pdf 22. Yarovsky, I.; Evans, E., (2002) Computer simulation of structure and properties of crosslinked polymers: application to epoxy resins. Polymer, 43(3): p. 963 23. Valavala, P.K., et al., (2007) Nonlinear multiscale modeling of polymer materials. International Journal of Solids and Structures. 44, (3-4), p. 1161 24. Mansoori, G. A.; George, T. F.; Zhang, G.; (2003) 'Measurement, Simulation and Prediction of Intermolecular Interactions and Structural Characterization of Organic Nanostructures' in Proceed.

of Conference on

Nanodevices and Systems, Nanotech, San Francisco, CA, p. 23 25. Ruberto, R.; Abramo, M. C.; Caccamo, C.; (2004) Molecular dynamics simulation study of binary fullerene mixtures Phys.

Rev.

B, 70 (15),

p.155413 26. Yamazaki, T.; (2011) Parametrization of coarse grained force elds for dynamic property of ethylene glycol oligomers/water binary mixtures arXiv:1112.4591v1, math_ph

27. Marrink, S.J.; Fuhrmans, M.; Risselada, H.J.; Periole, X.; (2008) The MARTINI force eld. In "Coarse graining of condensed phase and biomolecular systems", G. Voth ed., CRC press, Chapter 2. 28. Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S. J.; (2008) The MARTINI Coarse-Grained Force Field: Extension to Proteins Journal of Chemical Theory and Computation, 4 (5), p. 819 29. Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E.; (2008) GROMACS 4: Algorithms for Highly Ecient, Load-Balanced, and Scalable Molecular Simulation Journal of Chemical Theory and Computation, 4 (3), p. 435. 30. Akkermans, R.L.C.; Briels, W.J.; (2000) 'Coarse-grained dynamics of one chain in a polymer melt' J.Chem.Phys. 113 (15), p. 6409 31. Wolf, B. A.; (2003) Polymer-polymer interactions: consistent modeling in terms of chain connectivity and conformal response Macromolecular chemistry and physics, 207, p. 65 32. Lacombe, R. H.; Sanchez, I. C.; (1976) Statistical thermodynamics of uid mixtures The Journal of Physical Chemistry, 80 (23), 2568 33. Lacombe, R. H.; Sanchez, I. C.; (1978) Statistical Thermodynamics of Polymer Solutions Macromolecules, 11 (6), 1145 34. Nwabunma, D.; Kyu, T.; (1999) Phase Behavior of Mixtures of Low Molar Mass Nematic Liquid Crystal and in Situ Photo-Cross-Linked Polymer Network Macromolecules, 32 (3), p. 664 35. Maier, W.; Saupe, A.; Z. Naturforsch, (1959) A 14a, p. 882 36. Lemkul, J. A.; Allen, W. J.; Bevan, D. R.; (2010) Practical Considerations for Building GROMOS-Compatible Small-Molecule Topologies Journal of Chemical Information and Modeling, 50 (12), p. 2221 37. Fan, C.F.; Olafson, B.D.; Blanco, M.; Hsu, S.L.; (1992) Application of molecular simulation to derive phase diagrams of binary mixtures Macromolecules, 25 (14), p. 366 38. Huang, D. M.; Moule, A. J.; Faller, R.; (2011) Characterization of polymerfullerene mixtures for organic photovoltaics by systematically coarse-grained molecular simulations Fluid Phase Equilibria, 302(1-2), p. 21 39. Zhang, W. H.; Fu, B. X.; Seo, Y.; et al. (2002) Eect of methyl methacrylate/polyhedral oligomeric silsesquioxane random copolymers in compatibilization of polystyrene and poly(methyl methacrylate) blends Macromolecules, 35 (21), p. 8029 40. Ao, Z. M.; Jiang, Q.; (2006) Size eects on miscibility and glass transition temperature of binary polymer blend lms Langmuir, 22 (3), p. 1241 41. Icoz, D. Z.; Kokini, J. L.; (2007) Examination of the validity of the in terms of miscibility Flory-Huggins solution theory in dextran systems Carbohydrate Polymers, 68 (1), p. 59

42. Huynh, L.; Grant, J.; Leroux, J-C.; Delmas, P.; Allen, C.; (2007) Predicting the Solubility of the Anti-Cancer Agent Docetaxel in Small Molecule Excipients using Computational Methods Pharmaceutical Research, 25 (1), p. 147 43. Huang, D. M.; Faller, R.; Do, K.; Moulé, A. J.; (2010) Coarse-Grained Computer Simulations of Polymer/Fullerene Bulk Heterojunctions for Organic Photovoltaic Applications Journal of Chemical Theory and Computation, 6 (2), p. 526 44. Yang, C.; Hu, J. G.; Heeger, A. J.; (2006) Molecular Structure and Dynamics at the Interfaces within Bulk Heterojunction Materials for Solar Cells Journal of the American Chemical Society, 128 (36), p. 12007 45. Dantanarayana, V.; Huang, D. M.; Staton, J. A.; Moulé, A. J.; Faller, R.; (2012) Multi-Scale Modeling of Bulk Heterojunctions for Organic Photovoltaic Applications Third Generation Photovoltaics, Chapter 2, ISBN: 978-953-51-0304-2 46. Yang, C.; Zhou, E.; Miyanishi, S.; Hashimoto, K.; Tajima, K.; (2011) Preparation of Active Layers in Polymer Solar Cells by Aerosol Jet Printing ACS Applied Materials & Interfaces, 3 (10), p. 4053 47. Sumpter, B. G.; Meunier, V.; (2012) Can computational approaches aid in untangling the inherent complexity of practical organic photovoltaic systems?

J. Polym.

Sci.

B Polym.

Phys., 50 (11) p.

n/a (published

on-line) 48. Sekkal, W.; Aourag, H.; Certier, M.; (1999) The miscibility of C60xC70(1x) alloy at high temperature Physics Letters A, Atomic and Solid State Physics, 251 (2), p. 132 49. Girifalco, L. A.; (1992) Molecular properties of fullerene in the gas and solid phases The Journal of Physical Chemistry 96 (2), 858 50. Bartlett, P.; (1990) A model for the freezing of binary colloidal hard spheres J. Phys.: Condens. Matter 2, p. 4979 51. Kniaz, K.; Fischer, J. E.; McGhie, A. R.; Girifalco, L. A.; Strongin, R. M.; Smith, A. B.; (1995) Fullerene alloys Solid State Comm., 96, p. 739 52. Hume-Rothery, W.; Smallman, R. E.; Haworth, C.W.; The structure of metals and alloys (1969) Metals & Metallurgy Trust