The MesoDyn project: software for mesoscale chemical engineering

Journal of Molecular Structure (Theochem) 463 (1999) 139–143

The MesoDyn project: software for mesoscale chemical engineering P. Altevogt a,*, O.A. Evers b, J.G.E.M. Fraaije, N.M. Maurits, B.A.C. van Vlimmeren c a

IBM Science Heidelberg, Vangerowstr. 18, 69115 Heidelberg, Germany b BASF, 67056 Ludwigshafen, Germany c University of Groningen, Department of Biphysical Chemistry, Nijenborgh 4, 9747 AG Groningen, Netherlands

Abstract We describe a new class of phenomenological mesoscopic models to simulate the phase separation dynamics in three dimensional complex liquids, based on dynamic density functional methods. These models are generalizations of time-dependent Ginzburg–Landau models and contain a molecular description of the liquids in the free energy functional. Possible applications are in process industries (HIPS, paints, detergents, surfactants,…), petroleum industries (oil recovery), pharmaceutical industries (drug delivery) and consumer product industries (food processing, cosmetics). We show the results of the simulation of the microphase behaviour of aqueous PL64 solutions. The work described here was mainly done within ESPRIT research projects funded by the European Union. 䉷 1999 Published by Elsevier Science B.V. All rights reserved. Keywords: Mesoscale structures; Modeling; Software

1. Introduction Mesoscale 1 structures are of utmost importance during the production processes of many materials, e.g. polymer blends, block-copolymer systems, surfactant aggregates in detergent materials, and latex particles and drug delivery systems. Although there exists a wealth of experimental data about these structures, computer simulations based on mathematical models providing a deeper understanding, especially of their time development, are still rare. Mesoscopic dynamics models receive increasing attention as they form a bridge between fast molecular kinetics and slow thermodynamic relaxation of macroscale properties [1–5]. In line with these devel* Corresponding author. 1 Mesoscale stands for ca. 10–100 nm.

opments, the European Union-funded ESPRIT project, MesoDyn, aims at developing and implementing computational methods for the rational analysis and design of mesoscale structures in polymer liquids. The MesoDyn project is a collaboration between BASF (Germany), IBM Deutschland Informations systeme GmbH (Germany), Molecular Simulations Limited (UK), Shell Chemicals BV (Netherlands), Norsk Hydro A.S.A. (Norway) and the University of Groningen (Netherlands). The starting point for the mesoscale simulations is a coarse grained model for the diffusive and hydrodynamic phenomena in phase-separation dynamics. The thermodynamic forces are obtained via a mean-field density functional theory, assuming a Gaussian chain as a molecular model. The melt dynamics are described by a set of stochastic partial differential equations (functional Langevin equations) for polymer diffusion. Noise sources, with correlations

0166-1280/99/$ - see front matter 䉷 1999 Published by Elsevier Science B.V. All rights reserved. PII: S0166-128 0(98)00403-5

140

P. Altevogt et al. / Journal of Molecular Structure (Theochem) 463 (1999) 139–143

dictated by the fluctuation-dissipation theorem, introduce the thermal fluctuations. The numerical calculations involve the time-integration of functional Langevin equations, given an implicit Gaussian density functional expression for the intrinsic chemical potentials. For industrial-relevant polymer systems, these calculations are highly resource consuming (CPU, memory and bandwidth between CPU and memory). This makes the use of high performance computing (HPC) tools, especially multiprocessor systems with distributed memory, almost mandatory. The problems in developing the MesoDyn software can be split up into two categories. The first category consists of modeling problems, i.e. to find a mathematical model that represents the most important physical features of the mesoscale structures. The second category consists of algorithmic and implementation related issues, i.e. given a certain model, find the algorithms that solve the discretized problem accurately and efficiently. In this paper we present the current status of the model and the algorithms.

2. The model Dynamic density functional models consist of a free energy functional describing the physical systems under consideration and a set of stochastic partial differential equations determining their time evolution. A prototype of such a model [5,6] will be described in this section.

2.1. The free energy The model used in the MesoDyn project consists of beads of various types I,J,… with interactions described by harmonic oscillator potentials for the intra molecular interactions (Gaussian chain) [7] and a mean field potential for all other interactions. The goal of a simulation is to predict the time evolution of the density distributions for the bead types I,J,…, i.e. to describe the morphologies of the liquid. The ideal free energy functional F id describing n non-interacting Gaussian chains consisting of N

beads is given by [6]: F id ‰rŠ ˆ ⫺b⫺1 nlnF ⫹ b⫺1 lnn! ⫺

XZ

UI …r†rI …r†dr

I

…1†

FˆC

⫺1

Z

⫺b‰H G ⫹

e

N P

Us …Rs †Š

sˆ1

dR1 …dRN

…2†

with external potentials, UI, HG ˆ

N 3b⫺1 X …Rs ⫺ Rs⫺1 †2 2a2 sˆ2

CˆL

3

2pa2 3

…3†

! 3 …N⫺1† 2

…4†

s h2 b Lˆ 2pm

…5†

bead mass, m and the Planck constant, h. The density distribution in the system is given by: Z rI ‰UŠ…r† ˆ C 0 d…r

⫺ RI †e

⫺b‰H G …R1 ;…;RN †⫹

N P

Us …Rs †Š

sˆ1

dR1 …dRN …6†

0

with a normalisation constant, C ; for details see Ref. [6]. The mapping between {U} and {r} is bijective [8]. All other interactions between the beads are modeled by the mean field free energy: 1 X ZZ eIJ …兩r ⫺ r 0 兩†rI …r†rJ …r 0 †drdr 0 …7† F mf ‰rŠ ˆ 2 I;J with the Gaussian kernel:  3 2 ⫺ 3 2 …r⫺r 0 †2 3 eIJ …兩r ⫺ r 0 兩† ˆ eoIJ e 2a 2 2p a

…8†

The total free energy is defined as: F‰rŠ ˆ F id ‰rŠ ⫹ F mf ‰rŠ:

…9†

2.2. Dynamics The dynamics of the model are described by the

P. Altevogt et al. / Journal of Molecular Structure (Theochem) 463 (1999) 139–143

following set of diffusion equations (generalized timedependent Ginzburg–Landau model [9]): N Z X 2rI …r; t† dF 0 …r ; t†dr 0 ⫹ hI …r; t† ˆ DIJ …r; r 0 ; t† 2t dr J Jˆ1

…10† with diffusion operator, DIJ and noise, hI with correlations: 具hI …r; t† ˆ 0典

…11†

具hI …r; t†hJ …r 0 ; t 0 †典 ˆ ⫺2b⫺1 DIJ …r; r 0 ; t†d…t ⫺ t 0 †

…12†

The set of diffusion equations is closed by the expression for the free energy:

dF dF mf ˆ ⫺UI ⫹ drI drI

…13†

3. Numerics and implementation For the implementation of MesoDyn, we need a solver for the system of stochastic differential equations, the equations specifying the noise and the closure relation Eqs. (10)–(13) 2. In this section, we present an outline of our algorithm and its implementation using a finite difference scheme on a simple cubic grid [6,13]. 3.1. The algorithm Input: density fields {rI} and intrinsic chemical potentials {dF=drI }. Calculate {dF mf =drI } using {UI ˆ ⫺…dF=drI † ⫹ dF mf =drI }. Crank–Nicholson scheme:

…14†

with diffusion constant D. Non-local coupling has also been studied since it describes the physics of many systems more appropriately (see [6] and references therein). 2.3. Remarks The Gaussian chain model, time-dependent Ginzburg–Landau theories and density functional methods as described by Eqs. (9)–(13) are well established, but their combination is one model (MesoDyn) and the implementation of this model in three dimensions is unique. An important distinguishing feature of MesoDyn compared with classic Ginsburg–Landau models [9, 10] is the introduction of a molecular model. This allows for much more complex density distribution patterns (morphologies) to be generated during the time integration of the model. This is essential, e.g. in the simulation of microphase separation in polymer melts [6,11,12]. MesoDyn is a prototype for a whole new class of mesoscopic models. Studies to incorporate non-local diffusion operators, enhanced molecular descriptions, electrostatic effects, hydrodynamics and chemical reactions are in progress.

Eq.

(7)

and

1. Solve non-linear equations in {rI}, {dF=drI } resulting from discretisation of the differential equations. 2. Check for convergence and exit or go to step 3. 3. Select new {UI} using steepest descent or conjugate gradient methods. 4. Calculate {rI} [Eq. (6)] by applying an integration scheme based on Green propagators, {dF mf =drI } [Eq. (7)] and {dF=drI } [Eq. (13)] and go back to step 1.

If a local coupling approximation is used, the diffusion operator is given by: DIJ …r; r 0 † ˆ bDdIJ d…r ⫺ r 0 †7r ·rJ …r†7r 0

141

3.2. Implementation issues Implementation of steps 1 and 4 of the abovedescribed algorithm is accomplished by so called stencil operations [13]. These operations consist of loading fields (of densities or potentials) from memory into the CPU, calculating at each grid point a linear combination of these fields defined at neighbouring grid points and writing the resulting fields back to memory. The stencil operations are the numerical intensive kernels of MesoDyn and their efficient implementation on computer systems that should provide a good floating point performance, as well as a high bandwidth between memory and CPU, 2

To our knowledge, almost no mathematical theory exists for the numerics of this combination of stochastic partial differential and integral equations.

142

P. Altevogt et al. / Journal of Molecular Structure (Theochem) 463 (1999) 139–143

Fig. 3. Propylene oxide (PO) isosurface representations of aqueous PL64 solutions. The isosurfaces represent volume fractions of PO (uPO): (a) 70% PL64, isolevel uPO ˆ 0.59; (b) 60% PL64, isolevel uPO ˆ 0.50; (c) 55% PL64, isolevel PL64, isolvel uPO ˆ 0.46; and (d) 50% PL64, isolevel uPO ˆ 0.42. The colors denote the volume fraction of the ethylene oxide beads at the surface.

is crucial for the overall performance of a simulation. Considering, moreover, the demands of industrial applications (e.g. it should be possible to simulate the time development of an industrial relevant polymer blend for at least 1 msec within a few days) and the computer resources required for such simulations 3, the implementation of MesoDyn has been designed for multiprocessor systems of the MIMD type with distributed memory [14]. The parallelization of MesoDyn for such computer systems has been done by domain decomposition methods using the message passing programming model, as implemented in the MPI library [15]. The single node performance of the numerical kernels (especially of the stencil operations) has been optimized for RISC architectures with a memory 3 A simulation of a system consisting of 100 bead types on a 128 3 lattice would require more than 4.0 GB of main memory and would last many weeks on a high end workstation like, e.g. IBM RS/6000 590.

hierarchy by loop unrolling, blocking [14] and by exploiting the isotropy of the stencil operations to minimize the data traffic between processor and memory [13].

4. Applications Possible applications of MesoDyn include the simulation of the time-evolution of micelle formation (Evers, pers. comm.) (e.g. in Latex-based products), drug delivery systems, cosmetic products, oil recovery, food processing, formation of mesoscale structures in high impact polysterenes (HIPS) and coatings. Recently, MesoDyn has been used to simulate the microphase behaviour of aqueous PL64 [(ethylene oxide)13 (propylene oxide)30 (ethylene oxide)13] solutions to obtain the phase diagram and the structure factors, including their time development [12]. Fig. 1 shows the PL64 morphologies at several volume fractions.

P. Altevogt et al. / Journal of Molecular Structure (Theochem) 463 (1999) 139–143

5. Conclusions MesoDyn is based on a new combination of wellestablished concepts of polymer science and soft condensed matter physics. Its implementation uses advanced methods from numerics and computer science. Due to the enormous amount of possible industrial-relevant applications of MesoDyn, it has the potential to open up new areas for computer simulations in industry. To our knowledge, in combining all of the above mentioned features, MesoDyn is world-wide unique. Important next steps in future developments of MesoDyn will be the inclusion of new features (e.g. chemical reactions,…), investigating its parametrization and further experimental validation.

[4] [5]

[6]

[7] [8]

[9] [10]

Acknowledgements The work described in this article has been supported by the European Union in the framework of the ESPRIT projects EP8328 (CAESAR) and EP22685 (MesoDyn). References [1] O.T. Valls, J.E. Farrell, Spinodal decomposition in a threedimensional fluid model, Phys. Rev. E 47 (1993) R36–R39. [2] L. Ramı´rez-Piscina, A. Hernande´z-Machado, J.M. Sancho, Numerical algorithm for Ginzburg–Landau equations with multiplicative noise: application to domain growth, Phys. Rev. B 48 (1993) 125–131. [3] T. Kawakatsu, K. Kawasaki, M. Furusaka, H. Okabayashi, T. Kanaya, Late stage dynamics of phase separation processes of

[11]

[12]

[13]

[14]

[15]

143

binary mixtures containing surfactants, J. Chem. Phys. 99 (10) (1993) 8200–8217. A. Shinozaki, Y. Oono, Spinodal decomposition in 3-space, Phys. Rev. E 48 (1993) 2622–2654. J.G.E.M. Fraaije, Dynamic density functional theory for micro-phase separation kinetics of block copolymer melts, J. Chem. Phys. 99 (1993) 9202–9212. J.G.E.M. Fraaije, B.A.C. van Vilmmeren, N.M. Maurits, M. Postma, O.A. Evers, C. Hoffmann, P. Altevogt, G. GoldbeckWood, The dynamic mean-field density functional method and its application to the mesoscopic dynamics of quenched block copolymer melts, J. Chem. Phys. 106 (10) (1997) 4260– 4269. M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Oxford Science Publications, Oxford, 1986. J.T. Chayes, L. Chayes, On the validity of the inverse conjecture in classical density functional theory, J. Stat. Phys. 36 (3/ 4) (1984) 471–488. P.M. Chakin, T.C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, Cambridge, 1995. M.C. Cross, P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65 (3) (1993) 51–1112. R. Hasegawa, M. Doi, Adsorption dynamics. Extension of self-consistent field theory to dynamics. Extension of selfconsistent field theory to dynamical problems, Macromolecules 30 (1997) 3086–3089. B.A.C. van Vlimmeren, N.M. Maurits, A.V. Zvelindovsky, G.J.A. Sevink, J.G.E.M. Fraaije. Micro-phase separation kinetics in concentrated aqueous solution of triblock polymer surfactant (EP)13 (PO)30 (EO)13: an application of dynamic mean-field density functional theory (submitted). N.M. Maurits, P. Altevogt, O.A. Evers, J.G.E.M. Fraaije, Simple numerical quadrature rules for Gaussian chain polymer density functional calculations in 3D and implementation on parallel platforms, Comp. Pol. Sci. 6 (1996) 1. J.L. Hennessy, D.A. Patterson, Computer Architecture: A Quantitative Approach. Morgan Kaufmann, San Francisco, 1996. Message Passing Interface Forum. MPI: A Message Passing Interface Standard, http://www.mcs.anl.gov/mpi/index.html, 1997.