Efficient and realistic simulation of phase coexistence.

We show how an existing concurrent multi-scale method named hybrid particle field-molecular dynamics (hPF-MD) can be adapted to enable the simulation of structure and/or structural dynamics in compressible systems. Implementing such new equations of state (EOS) into hPF-MD, while conserving the efficiency associated with treating intermolecular interactions in a continuum fashion, opens this method up to describe a new class of phenomena in which non-uniform densities play a role, for example, evaporation and crystallization. We carefully consider how compressible hPF-MD compares to its mean-field counterpart for two particular EOS, adopted from the Cell Model for polymers and the Carnahan-Starling expression for hard spheres. Here, we performed a very basic analysis for a single-component system, focusing on the significance of various particle-based parameters and the particle-to-field projection. Our results illustrate the key role of the particle density per field grid cell and show that projection based on a Gaussian kernel is preferred over the standard cloud-in-cell projection. They also suggest that the behavior of hPF-MD close to the critical point is non-classical, i.e., in agreement with a critical exponent for a pure particle description, despite the mean-field origin of the method.

Combining cell-based hydrodynamics with hybrid particle-field simulations: efficient and realistic simulation of structuring dynamics.

We have extended an existing hybrid MD-SCF simulation technique that employs a coarsening step to enhance the computational efficiency of evaluating non-bonded particle interactions. This technique is conceptually equivalent to the single chain in mean-field (SCMF) method in polymer physics, in the sense that non-bonded interactions are derived from the non-ideal chemical potential in self-consistent field (SCF) theory, after a particle-to-field projection. In contrast to SCMF, however, MD-SCF evolves particle coordinates by the usual Newton's equation of motion. Since collisions are seriously affected by the softening of non-bonded interactions that originates from their evaluation at the coarser continuum level, we have devised a way to reinsert the effect of collisions on the structural evolution. Merging MD-SCF with multi-particle collision dynamics (MPCD), we mimic particle collisions at the level of computational cells and at the same time properly account for the momentum transfer that is important for a realistic system evolution. The resulting hybrid MD-SCF/MPCD method was validated for a particular coarse-grained model of phospholipids in aqueous solution, against reference full-particle simulations and the original MD-SCF model. We additionally implemented and tested an alternative and more isotropic finite difference gradient. Our results show that efficiency is improved by merging MD-SCF with MPCD, as properly accounting for hydrodynamic interactions considerably speeds up the phase separation dynamics, with negligible additional computational costs compared to efficient MD-SCF. This new method enables realistic simulations of large-scale systems that are needed to investigate the applications of self-assembled structures of lipids in nanotechnologies.

Erratum. Computational insight in the role of fusogenic lipopeptides at the onset of liposome fusion.

We performed an extensive computational study to obtain insight in the molecular mechanisms that take place prior to membrane fusion. We focused on membrane-anchored hybrid macromolecules (lipid­polymer­oligopeptide) that mimic biological SNARE proteins in terms of liposome fusion characteristics [H. Robson Marsden et al., 2009]; efficient micro-second simulation was enabled by combining validated MARTINI force fields for the molecular building blocks in coarse-grained molecular dynamics (CGMD). We find that individual peptide domains in the hybrid macromolecules bind and partially integrate parallel to the membrane surface, in agreement with experimental findings. By varying several experimental design parameters, we observe that peptide domains remain in the solvent phase only in two cases: (1) for solitary lipopeptides (low concentration), below a threshold area per lipid in the membrane, and (2) when the lipopeptide concentration is high enough for the peptide domains to self-assemble into tetrameric homo-complexes. The peptide-membrane binding is not affected by solvent-induced peptide unfolding, which we mimicked by relaxing the usual MARTINI helix constraints. Remarkably, in this case, a reverse transition to a helical secondary structure is observed after binding, highlighting the role of the membrane as a template (partitioning-folding coupling). Our findings undermine the current view of the initial stages towards fusion, in which membranes are thought to be kept in close apposition via dimerization of individual complementary peptides in the solvent phase. Although we did not study actual fusion, our simulations show that the formation of homomers, which is suppressed in experimental peptide pair design and therefore believed to be insignificant for fusion, by peptides anchored to the same membrane does play a key role in this locking mechanism and potentially also in membrane destabilization that precedes fusion.

Rigorous embedding of cell dynamics simulations in the Cahn-Hilliard-Cook framework: Imposing stability and isotropy.

We have rigorously analyzed the stability of the efficient cell dynamics simulations (CDS) method by making use of the special properties of the local averaging operator 〈〈*〉〉-* in matrix form. Besides resolving a theoretical issue that has puzzled many over the past three decades, this analysis has considerable practical value: It relates CDS directly to finite-difference approximations of the Cahn-Hilliard-Cook equations and provides a straightforward recipe for replacing the original two- or three-dimensional (2D or 3D) averaging operators in CDS by an equivalent (in terms of stability) discrete Laplacian with superior isotropy and scaling behavior. As such, we open up a route to suppress the unphysical reflection of the computational grid in CDS results (grid artifacts). We found that proper rescaling of discrete Laplacians, needed to employ them in CDS, is equivalent to introducing a well-chosen time step in CDS. In turn, our analysis provides stability conditions for phase-field simulations based on the Cahn-Hilliard-Cook equations. Subsequently, we have quantitatively compared the isotropy and scaling behavior of several discrete 2D or 3D Laplacians, thereby extending the significance of this work to general field-based methodology. We found that all considered discrete Laplacians have equivalent scaling behavior along the Cartesian directions. In addition, and somewhat surprisingly, known "isotropic" discrete Laplacians, i.e., isotropic up to fourth order in |k|, become quite anisotropic for larger wave vectors, whereas "less isotropic" discrete Laplacians (second order) are only slightly anisotropic on the whole |k| range. We identified a hard limit to the accuracy with which the discrete Laplacian can emulate the two important properties of the optimal (continuum) Laplacian, as an improvement of the isotropy, by introducing additional points to the stencil, will negatively affect the scaling behavior. Within this limitation, the discrete compact Laplacians in the DnQm class known from lattice hydrodynamics, D2Q9 in 2D and D3Q19 in 3D, are found to be optimal in terms of isotropy. However, by being only slightly anisotropic on the whole range and enabling larger time steps, the discrete Laplacians that relate to the local averaging operator of Oono and Puri (2D) and Shinozaki and Oono (3D) as well as the less familiar 3D discrete BvV Laplacian developed for dynamic density functional theory are valid alternatives.

Computational insight in the role of fusogenic lipopeptides at the onset of liposome fusion.

We performed an extensive computational study to obtain insight in the molecular mechanisms that take place prior to membrane fusion. We focused on membrane-anchored hybrid macromolecules (lipid-polymer-oligopeptide) that mimic biological SNARE proteins in terms of liposome fusion characteristics [H. Robson Marsden et al., 2009]; efficient micro-second simulation was enabled by combining validated MARTINI force fields for the molecular building blocks in coarse-grained molecular dynamics (CGMD). We find that individual peptide domains in the hybrid macromolecules bind and partially integrate parallel to the membrane surface, in agreement with experimental findings. By varying several experimental design parameters, we observe that peptide domains remain in the solvent phase only in two cases: (1) for solitary lipopeptides (low concentration), below a threshold area per lipid in the membrane, and (2) when the lipopeptide concentration is high enough for the peptide domains to self-assemble into tetrameric homo-complexes. The peptide-membrane binding is not affected by solvent-induced peptide unfolding, which we mimicked by relaxing the usual MARTINI helix constraints. Remarkably, in this case, a reverse transition to a helical secondary structure is observed after binding, highlighting the role of the membrane as a template (partitioning-folding coupling). Our findings undermine the current view of the initial stages towards fusion, in which membranes are thought to be kept in close apposition via dimerization of individual complementary peptides in the solvent phase. Although we did not study actual fusion, our simulations show that the formation of homomers, which is suppressed in experimental peptide-pair design and therefore believed to be insignificant for fusion, by peptides anchored to the same membrane does play a key role in this locking mechanism and potentially also in membrane destabilization that precede fusion.

Efficient solvent-free dissipative particle dynamics for lipid bilayers.

We rigorously derived effective potentials for solvent-free DPD simulation of lipid bilayers. The derivation relies on an earlier developed hybrid particle/field method and is based on the idea that the solvent is always in local equilibrium on a coarse-grained time scale, given the instantaneous templates set by the self-assembly structure. By relating the parameters in the effective implicit-solvent potentials directly to the lipid-solvent interactions and membrane properties for the explicit solvent DPD model, we constitute an efficient and general procedure for reformulating any DPD membrane model in an implicit-solvent form. Here, we determined these membrane properties for two existing DPD models, via an analysis of membrane fluctuation spectra. Equivalent single-processor implicit- and explicit-solvent calculations show the trade-mark of implicit solvent simulation: a 20-fold reduction of the total simulation time for a system containing 92% solvent. This increased efficiency enabled us to realistically simulate the spontaneous formation of a ∼20 nm diameter vesicle on a single processor overnight. We believe that this work will contribute to an enhanced computational study of large vesicles and thus a better understanding of experimental liposome dynamics.

Stochastic quasi-Newton method: application to minimal model for proteins.

Knowledge of protein folding pathways and inherent structures is of utmost importance for our understanding of biological function, including the rational design of drugs and future treatments against protein misfolds. Computational approaches have now reached the stage where they can assess folding properties and provide data that is complementary to or even inaccessible by experimental imaging techniques. Minimal models of proteins, which make possible the simulation of protein folding dynamics by (systematic) coarse graining, have provided understanding in terms of descriptors for folding, folding kinetics, and folded states. Here we focus on the efficiency of equilibration on the coarse-grained level. In particular, we applied a new regularized stochastic quasi-Newton (S-QN) method, developed for accelerated configurational space sampling while maintaining thermodynamic consistency, to analyze the folding pathway and inherent structures of a selected protein, where regularization was introduced to improve stability. The adaptive compound mobility matrix B in S-QN, determined by a factorized secant update, gives rise to an automated scaling of all modes in the protein, in particular an acceleration of protein domain dynamics or principal modes and a slowing down of fast modes or "soft" bond constraints, similar to lincs/shake algorithms, when compared to conventional Langevin dynamics. We used and analyzed a two-step strategy. Owing to the enhanced sampling properties of S-QN and increased barrier crossing at high temperatures (in reduced units), a hierarchy of inherent protein structures is first efficiently determined by applying S-QN for a single initial structure and T=1>T(θ), where T(θ) is the collapse temperature. Second, S-QN simulations for several initial structures at very low temperature (T=0.01

Stochastic quasi-Newton molecular simulations.

We report a new and efficient factorized algorithm for the determination of the adaptive compound mobility matrix B in a stochastic quasi-Newton method (S-QN) that does not require additional potential evaluations. For one-dimensional and two-dimensional test systems, we previously showed that S-QN gives rise to efficient configurational space sampling with good thermodynamic consistency [C. D. Chau, G. J. A. Sevink, and J. G. E. M. Fraaije, J. Chem. Phys. 128, 244110 (2008)]. Potential applications of S-QN are quite ambitious, and include structure optimization, analysis of correlations and automated extraction of cooperative modes. However, the potential can only be fully exploited if the computational and memory requirements of the original algorithm are significantly reduced. In this paper, we consider a factorized mobility matrix B=JJ(T) and focus on the nontrivial fundamentals of an efficient algorithm for updating the noise multiplier J . The new algorithm requires O(n2) multiplications per time step instead of the O(n3) multiplications in the original scheme due to Choleski decomposition. In a recursive form, the update scheme circumvents matrix storage and enables limited-memory implementation, in the spirit of the well-known limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method, allowing for a further reduction of the computational effort to O(n). We analyze in detail the performance of the factorized (FSU) and limited-memory (L-FSU) algorithms in terms of convergence and (multiscale) sampling, for an elementary but relevant system that involves multiple time and length scales. Finally, we use this analysis to formulate conditions for the simulation of the complex high-dimensional potential energy landscapes of interest.

Improved configuration space sampling: Langevin dynamics with alternative mobility.

We present a new and efficient method for determining optimal configurations of a large number (N) of interacting particles. We use a coarse-grained stochastic Langevin equation in the overdamped limit to describe the dynamics of this system and replace the standard mobility by an effective space dependent inverse Hessian correlation matrix. Due to the analogy of the drift term in the Langevin equation and the update scheme in Newton's method, we expect accelerated dynamics or improved convergence in the convex part of the potential energy surface Phi. The stochastic noise term, however, is not only essential for proper thermodynamic sampling but also allows the system to access transition states in the concave parts of Phi. We employ a Broyden-Fletcher-Goldfarb-Shannon method for updating the local mobility matrix. Quantitative analysis for one and two dimensional systems shows that the new method is indeed more efficient than standard methods with constant effective friction. Due to the construction, our effective mobility adapts high values/low friction in configurations which are less optimal and low values/high friction in configurations that are more optimal.

Block copolymers confined in a nanopore: pathfinding in a curving and frustrating flatland.

We have studied structure formation in a confined block copolymer melt by means of dynamic density functional theory. The confinement is two dimensional, and the confined geometry is that of a cylindrical nanopore. Although the results of this study are general, our coarse-grained molecular model is inspired by an experimental lamella-forming polysterene-polybutadiene diblock copolymer system [K. Shin et al., Science 306, 76 (2004)], in which an exotic toroidal structure was observed upon confinement in alumina nanopores. Our computational study shows that a zoo of exotic structures can be formed, although the majority, including the catenoid, helix, and double helix that were also found in Monte Carlo nanopore studies, are metastable states. We introduce a general classification scheme and consider the role of kinetics and elongational pressure on stability and formation pathway of both equilibrium and metastable structures in detail. We find that helicity and threefold connections mediate structural transitions on a larger scale. Moreover, by matching the remaining parameter in our mesoscopic method, the Flory-Huggins parameter chi, to the experimental system, we obtain a structure that resembles the experimental toroidal structure in great detail. Here, the most important factor seems to be the roughness of the pore, i.e., small variations of the pore radius on a scale that is larger than the characteristic size in the system.

Scaling behavior of the reorientation kinetics of block copolymers exposed to electric fields.

We have followed the reorientation kinetics of various block copolymer solutions exposed to an external electric DC field. The characteristic time constants follow a power law indicating that the reorientation is driven by a decrease in electrostatic energy. Moreover, the observed exponent suggests an activated process in line with the expectations for a nucleation and growth process. When properly scaled, the data collapse onto a single master curve spanning several orders of magnitude both in reduced time and in reduced energy. The power law dependence of the rate of reorientation derived from computer simulations based on dynamic density functional theory agrees well with the experimental observations. First experiments in AC electric fields at sufficiently high frequencies confirm the notion that the reorientation process is dominated by differences in the dielectric constants rather than by mobile ions.

Dynamics of terrace formation in a nanostructured thin block copolymer film.

We have used dynamic self-consistent field (DSCF) theory to investigate the structural evolution of an ABA block copolymer thin film placed between a solid substrate and a free surface. In line with the few existing theoretical studies for pure homopolymers and mixtures, the free interface is introduced by a void component. In our calculations, the free surface experiences surface roughening and eventually the formation of terraces, as in the experiments. The kinetic pathway of the microstructures was compared to findings of an existing detailed experimental study (Knoll, A.; Lyakhova, K. S.; Horvat, A.; Krausch, G.; Sevink, G. J. A.; Zvelindovsky, A. V.; Magerle, R. Nat. Mater. 2004, 3, 886) and was found to be equivalent in detail. This corroborates our assumption in this earlier work that the pathway due to changing film thickness is similar to a pathway due to changing surface energetics. Moreover, our calculations show for the first time that microstructural transitions are a driving force of polymer/air interface curving and the formation of terraces.

Influence of initial order on the microscopic mechanism of electric field induced alignment of block copolymer microdomains.

We investigate the mechanism of microdomain orientation in concentrated block copolymer solutions exposed to a dc electric field by in situ synchrotron small-angle X-ray scattering (SAXS). As a model system, we use concentrated solutions of a lamellar polystyrene-b-polyisoprene block copolymer in toluene. We find that both the microscopic mechanism of reorientation and the kinetics of the process strongly depend on the initial degree of order in the system. In a highly ordered lamellar system with the lamellae being aligned perpendicular to the electric field vector, only nucleation and growth of domains is possible as a pathway to reorientation and the process proceeds rather slowly. In less ordered samples, grain rotation becomes possible as an alternative pathway, and the process proceeds considerably faster. The interpretation of our finding is strongly corroborated by dynamic self-consistent field simulations.

Orthogonal fields: a path to long-range three-dimensional order in block copolymers.

Large-scale computer simulations show that two orthogonal external fields can control the orientation of lamellar microdomains in diblock copolymers in three dimensions and lead to an enhanced long-range ordering.

Design of chimaeric polymersomes.

We discuss the development of hierarchical polymer particles, or variegated polymersome composites, in which at least two different components are phase separated within one polymersome chimaera. We briefly discuss the present status in experimental polymersome research, and then discuss a speculative design strategy, based on mesoscopic simulations with a dynamical variant of polymer self-consistent field theory (Mesodyn). The main conclusion is that the counter-intuitive co-assembly of demixing block copolymers is the key in controlling hierarchical structures on a mesoscopic scale. This is the classical paradox of a chimaera: the constituents live in the same scaffold, but apart. Block copolymers beyond a certain length will always split the assembly, and without further precautions, polymer based chimaerae are intrinsically unstable. To this end, we propose the application of a branched block copolymer as composite compatibilizer, glueing the separate domains together, and thereby stabilizing the chimaeric polymersome.

Direct imaging and mesoscale modelling of phase transitions in a nanostructured fluid.

The kinetics of phase transitions is essential for understanding pattern formation in structured fluids. These fluids play a key role in the morphogenesis of biological cells, and they are very common in pharmaceutical products and plastic materials. Until now, it has not been possible to follow phase transitions in structured fluids experimentally in real time and with high spatial resolution. Previous work has relied on static images and indirect experimental evidence from spatially averaging scattering experiments. Simulating the processes with computer models is a further challenge because of the multiple time and length scales involved. Our movies based on in situ scanning force microscopy show the time sequence of the elementary steps of a phase transition in a fluid film of block copolymer from the cylinder to the perforated lamella phase. The movies validate a versatile simulation model that gives physical insight into the nature of the process. Our approach provides a means of improving the study and understanding of pattern formation processes in nanostructured fluids. We expect a significant impact on nanotechnology where block copolymers serve as self-organized templates for the synthesis of inorganic nanostructured materials.

Kinetic pathways of sheared block copolymer systems derived from Minkowski functionals.

We employ Minkowski functionals to analyze the kinetics of pattern formation under an applied external shear flow. The considered pattern formation model describes the dynamics of phase separating block copolymer systems. For our purpose, we have chosen two block copolymer systems (a melt and a solution) that exhibit a hexagonal cylindrical morphology as an equilibrium structure. Our main objective is the determination of efficient choices for the treshold values that are required for the calculation of the Minkowski functionals. We find that a minimal set of two treshold values (one from which should be equal to an average density value and another to a higher density value) is sufficient to unraffle the phase separation kinetics. Given these choices, we focus on the influence of the degree of phase separation, and the instance at which the shear is applied, on the kinetic pathways. We also found a remarkable similarity of the time evolution of Euler characteristic and the segregation parameter for the average density choice.

Phase behavior in thin films of cylinder-forming ABA block copolymers: mesoscale modeling.

The phase behavior of cylinder-forming ABA block copolymers in thin films is modeled in detail using dynamic density functional theory and compared with recent experiments on polystyrene-block-polybutadiene-block-polystyrene triblock copolymers. Deviations from the bulk structure, such as wetting layer, perforated lamella, and lamella, are identified as surface reconstructions. Their stability regions are determined by an interplay between surface fields and confinement effects. Our results give evidence for a general mechanism governing the phase behavior in thin films of modulated phases.

Role of dissimilar interfaces in thin films of cylinder-forming block copolymers.

We study the effect of dissimilar interfaces on the phase behavior of cylinder forming block copolymers in thin films by means of dynamic density-functional theory. In this article, we show that dissimilarity of the interfaces induces hybrid structures. These structures appear when the surface fields at the two interfaces stabilize different surface structures and/or reconstructions. We propose a general classification of hybrid structures and give an unifying description of phase behavior of cylinder forming block copolymer films. Our results are consistent with experimental observations.