Mesoscale simulation of polymer reaction equilibrium: combining dissipative particle dynamics with reaction ensemble Monte Carlo. I. Polydispersed polymer systems

We present a mesoscale simulation technique, called the reaction ensemble dissipative particle dynamics (RxDPD) method, for studying reaction equilibrium of polymer systems. The RxDPD method combines elements of dissipative particle dynamics (DPD) and reaction ensemble Monte Carlo (RxMC), allowing for the determination of both static and dynamical properties of a polymer system. The RxDPD method is demonstrated by considering several simple polydispersed homopolymer systems. RxDPD can be used to predict the polydispersity due to various effects, including solvents, additives, temperature, pressure, shear, and confinement. Extensions of the method to other polymer systems are straightforward, including grafted, cross-linked polymers, and block copolymers. To simulate polydispersity, the system contains full polymer chains and a single fractional polymer chain, i.e., a polymer chain with a single fractional DPD particle. The fractional particle is coupled to the system via a coupling parameter that varies between zero (no interaction between the fractional particle and the other particles in the system) and one (full interaction between the fractional particle and the other particles in the system). The time evolution of the system is governed by the DPD equations of motion, accompanied by changes in the coupling parameter. The coupling-parameter changes are either accepted with a probability derived from the grand canonical partition function or governed by an equation of motion derived from the extended Lagrangian. The coupling-parameter changes mimic forward and reverse reaction steps, as in RxMC simulations.     

The ensemble folding dynamics of EF-hand domain in parvalbumin from a Monte Carlo simulation

The helix-loop-helix (i.e., EF-hand) is the most common motif in the superfamily of Ca\(^{2+}\)-binding proteins. Parvalbumins, as the classical EF-hand proteins, are associated with the muscle relaxation/contration, the calcium buffering etc. The previous researches focus more of interest on the ion coupled/decoupled mechanism using molecular dynamics (MD) computations. We developed the novel approach instead of MD, in which the Landau free energy was introduced to describe the protein in term of the skeletal C\(_\alpha \) chain. The unfolding and folding processes were simulated by the Glauber algorithm under the repeated heating and cooling cycles. The high-quality crystal structure (2PVB) was as the trigger of non-equilibrium dynamics simulation for parvalbumin-\(\beta \) (Parv). The evolution of three dimensional structure during the folding was displayed for the residues 8–64 fragment in Parv. The phenomenon of helix nucleation was observed, that was, the 3\(_{10}\) helix at residues 35–37 and the front of \(\alpha \)-helix at residues 40–50 were firstly formed in the folding process. We also found two potential misfoldings which were caused by the distortion of the local structures at residues 35–37, 45–50.     

Monte Carlo adiabatic simulation of equilibrium reacting systems: The ammonia synthesis reaction

We present an application of the recently developed Monte Carlo method for simulations at fixed total enthalpy [W. R. Smith, M. Lísal, Phys. Rev. E 66 (2002) 01114-1–01114-3], combined with the reaction ensemble Monte Carlo method, for the direct prediction of equilibrium reactive adiabatic processes. For the industrially important ammonia synthesis reaction in an adiabatic plug-flow reactor, we perform direct simulations of the equilibrium reaction temperature and the composition of the exit stream as a function of the temperature and pressure of the inlet stream. The chemical species of the system are represented by all-atom potentials with interaction parameters taken from the literature. The accuracy of the molecular model is validated by comparing simulation results with experimental data. We also compare the simulation results with a macroscopic thermodynamic model based on the Soave–Redlich–Kwong equation of state. The simulation results for the reaction conversion show very good agreement with available experimental data over a wide range of temperatures and pressures, whereas the corresponding results from the macroscopic thermodynamic model slightly deteriorate with increasing pressure. Based on these comparisons, the predicted values of the reaction temperature and composition of the exit stream from the simulations are more accurate than the corresponding predicted values from the macroscopic thermodynamic model.     

Monte Carlo simulation of a polymer-analogous reaction in a polymer blend

The autocatalytic polymer-analogous reaction A → B in a blend composed of two contacting layers of compatible homopolymers A and B is studied by numerical simulation using the dynamic continuum Monte Carlo method. The evolution of the numerical density of units A and units initially belonged to the chains of homopolymer A is investigated in the course of the reaction and interdiffusion. Local characteristics of the distribution of the homopolymer with respect to its composition and blocks A and B with respect to their length are calculated at different times. The dispersions of the above distributions are appreciably higher than the corresponding dispersion of the Bernoullian copolymer of the same average composition, despite the random character of the reaction. This effect can be provided by changes in the composition of the blend on the scale of the reacting chain as well as by the diffusive mixing of the above chains. For the products of the polymer-analogous reaction, the broadening of the compositional distribution is predicted also by the theoretical model, which describes interdiffusion in the reacting system on scales that are markedly greater than the size of a polymer chain.     

Application of kinetic Monte Carlo method to equilibrium systems: vapour-liquid equilibria

Kinetic Monte Carlo (kMC) simulations were carried out to describe the vapour-liquid equilibria of argon at various temperatures. This paper aims to demonstrate the potential of the kMC technique in the analysis of equilibrium systems and its advantages over the traditional Monte Carlo method, which is based on the Metropolis algorithm. The key feature of the kMC is the absence of discarded trial moves of molecules, which ensures larger number of configurations that are collected for time averaging. Consequently, the kMC technique results in significantly fewer errors for the same number of Monte Carlo steps, especially when the fluid is rarefied. An additional advantage of the kMC is that the relative displacement probability of molecules is significantly larger in rarefied regions, which results in a more efficient sampling. This provides a more reliable determination of the vapour phase pressure and density in case of non-uniform density distributions, such as the vapour-liquid interface or a fluid adsorbed on an open surface. We performed kMC simulations in a canonical ensemble, with a liquid slab in the middle of the simulation box to model two vapour-liquid interfaces. A number of thermodynamic properties such as the pressure, density, heat of evaporation and the surface tension were reliably determined as time averages.     

Dense orientationally ordered states of a single semiflexible macromolecule: an expanded ensemble Monte Carlo simulation

Using a coarse-grained model we perform a Monte Carlo simulation of the state behavior of an individual semiflexible macromolecule. Chains consisting of N = 256 and 512 monomer units have been investigated. A recently proposed enhanced sampling Monte Carlo technique for the bond fluctuation model in an expanded ensemble in four-dimensional coordinate space was applied. The algorithm allows one to accelerate the sampling of statistically independent three-dimensional conformations in a dense globular state. We found that the temperature of the intraglobular liquid-solid transition decreases with increasing chain stiffness. We have investigated the possible intraglobular orientationally ordered (i.e., liquid-crystalline) structures and obtained a diagram of states for chains consisting of N = 256 monomer units. This diagram contains regions of stability of coil, two spherical globules (liquid and solid), and rod-like globule conformations. Transitions between the globular states are rounded first-order ones since the states of liquid, solid, and cylinder-like globules do have different internal symmetry.     

Phase coexistence in heterogeneous porous media: a new extension to Gibbs ensemble Monte Carlo simulation method

The effect of confinement on phase behavior of simple fluids is still an area of intensive research. In between experiment and theory, molecular simulation is a powerful tool to study the effect of confinement in realistic porous materials, containing some disorder. Previous simulation works aiming at establishing the phase diagram of a confined Lennard-Jones-type fluid, concentrated on simple pore geometries (slits or cylinders). The development of the Gibbs ensemble Monte Carlo technique by Panagiotopoulos [Mol. Phys. 61, 813 (1987)], greatly favored the study of such simple geometries for two reasons. First, the technique is very efficient to calculate the phase diagram, since each run (at a given temperature) converges directly to an equilibrium between a gaslike and a liquidlike phase. Second, due to volume exchange procedure between the two phases, at least one invariant direction of space is required for applicability of this method, which is the case for slits or cylinders. Generally, the introduction of some disorder in such simple pores breaks the initial invariance in one of the space directions and prevents to work in the Gibbs ensemble. The simulation techniques for such disordered systems are numerous (grand canonical Monte Carlo, molecular dynamics, histogram reweighting, N-P-T+test method, Gibbs-Duhem integration procedure, etc.). However, the Gibbs ensemble technique, which gives directly the coexistence between phases, was never generalized to such systems. In this work, we focus on two weakly disordered pores for which a modified Gibbs ensemble Monte Carlo technique can be applied. One of the pores is geometrically undulated, whereas the second is cylindrical but presents a chemical variation which gives rise to a modulation of the wall potential. In the first case almost no change in the phase diagram is observed, whereas in the second strong modifications are reported.     

Communication: equation of state of hard oblate ellipsoids by replica exchange Monte Carlo

We implemented the replica exchange Monte Carlo technique to produce the equation of state of hard 1:5 aspect-ratio oblate ellipsoids for a wide density range. For this purpose, we considered the analytical approximation of the overlap distance given by Bern and Pechukas and the exact numerical solution given by Perram and Wertheim. For both cases we capture the expected isotropic-nematic transition at low densities and a nematic-crystal transition at larger densities. For the exact case, these transitions occur at the volume fraction 0.341, and in the interval 0.584-0.605, respectively.     

Excess energy and equation of state of fluids with hard-core potential models from a second-order Monte Carlo perturbation theory

Monte Carlo simulations in the NVT ensemble of the reference hard-sphere fluid have been performed to obtain the “exact” first- and second-order terms in the inverse temperature expansion of the free energy of fluids with hard-core potentials. The results have been used to obtain parametrizations of the free energy of fluids with Sutherland potentials with variable range as well as for a fluid with a hard-core Lennard–Jones potential. The results for the excess energy and the equation of state are compared with simulation data available in the literature for these fluids.     

Gibbs ensemble Monte Carlo simulation of supercritical CO2 adsorption on NaA and NaX zeolites

Adsorption of supercritical carbon dioxide on two kinds of zeolites with identical chemical composition but different pore structure (NaA and NaX) was studied using the Gibbs ensemble Monte Carlo simulation. The model frameworks for the two zeolites with SiAl ratio being unity have been chosen as the solid structures in the simulation. The adsorption behaviors of supercritical CO2 on the NaA and NaX zeolites, based on the adsorption isotherms and isosteric heats of adsorption, were discussed in detail and were compared with the available experimental results. A good agreement between the simulated and experimental results is obtained for both the adsorbed amount and the bulk phase density. The intermediate configurational snapshots and the radial distribution functions between zeolite and adsorbed CO2 molecules were collected in order to investigate the preferable adsorption locations and the confined structure behavior of CO2. The structure behaviors of the adsorbed CO2 molecules show various performances, as compared with the bulk phase, due to the confined effect in the zeolite pores.     

Vapor-liquid equilibria simulation and an equation of state contribution for dipole-quadrupole interactions

A systematic investigation on vapor-liquid equilibria (VLEs) of dipolar and quadrupolar fluids is carried out by molecular simulation to develop a new Helmholtz energy contribution for equations of state (EOSs). Twelve two-center Lennard-Jones plus point dipole and point quadrupole model fluids (2CLJDQ) are studied for different reduced dipolar moments micro*2=6 or 12, reduced quadrupolar moments Q*2=2 or 4 and reduced elongations L*=0, 0.505, or 1. Temperatures cover a wide range from about 55% to 95% of the critical temperature of each fluid. The NpT+test particle method is used for the calculation of vapor pressure, saturated densities, and saturated enthalpies. Critical data and the acentric factor are obtained from fits to the simulation data. On the basis of this data, an EOS contribution for the dipole-quadrupole cross-interactions of nonspherical molecules is developed. The expression is based on a third-order perturbation theory, and the model constants are adjusted to the present 2CLJDQ simulation results. When applied to mixtures, the model is found to be in excellent agreement with results from simulation and experiment. The new EOS contribution is also compatible with segment-based EOS, such as the various forms of the statistical associating fluid theory EOS.     

Monte Carlo simulation of solute extraction via supercritical carbon dioxide from poly(ethylene glycol)

The partitioning of ethylbenzene between poly(ethylene glycol) (PEG) and supercritical carbon dioxide was studied at 308.15, 328.15 and 348.15 K and 10, 15.5 and 20 MPa with PEG-400, 600 and 900 using Monte Carlo molecular simulation. The effect of a cosolvent was also studied with either 5% ethane or 5% n-octane added. Ethylbenzene favored the supercritical phase most when the density was highest, and while ethane had little effect, the addition of n-octane increased the amount of solute dissolved in carbon dioxide. Increasing polymer molecular weight led to more solute in the PEG-rich phase. This coincides with a higher amount of dissolved carbon dioxide that preferentially solvates ethylbenzene.     

Monte Carlo simulation of the interchain exchange reaction in a blend of incompatible polymers

The interchain exchange reaction in a blend composed of two contacting layers of incompatible A and B homopolymers was simulated by means of the dynamic off-lattice Monte Carlo method. The evolution of local molecular-mass and block mass distributions, depending on the effective temperature and the reaction rate, was studied for the first time. It was shown that the components interpenetrate as the copolymer forms in the interphase layer and the average block length decreases below a certain, temperature-dependent value. The state of dynamic equilibrium, whose characteristics are determined mainly by temperature, is established in the system. The time of establishment of equilibrium and the intensity of compatibilization at the early steps of the process are controlled by the rate of the reaction. The results of the study allow the contribution of the reaction to the interchange processes to be evaluated.     

Kinetic analysis of the cross reaction between dithioester and alkoxyamine by a Monte Carlo simulation

A model reaction of dithioester and alkoxyamine is proposed to probe the reversible addition–fragmentation chain transfer (RAFT) process. The kinetics of the model reaction is analyzed and compared with that of pure alkoxyamine homolysis with a Monte Carlo simulation. Although the pure alkoxyamine obeys the law of persistent radical effect, the model reaction results in higher concentration of the persistent radical during the main period of the reaction. However, for a very fast RAFT process or a very low addition rate constant, the time dependence of the persistent radical concentration is quite close to that of pure alkoxyamine. Furthermore, the cross termination between the intermediate and alkyl radicals causes a retardation effect for the model reaction when the intermediate is relatively long‐lived. The Monte Carlo simulation indicates that it is feasible to measure the individual rate constants of the RAFT process, such as the rate constant of addition, with a large excess of alkoxyamine. In addition, the special feature of the system with different leaving groups in the alkoxyamine and dithioester is also discussed. © 2006 Wiley Periodicals, Inc. J Polym Sci Part A: Polym Chem 45: 374–387, 2007     

Monte Carlo simulation of polymer analogous reaction in confined conditions: effects of ordering

The influence of energetic parameters of the interchain homo- and heterocontacts on a local ordering of Bernoullian copolymers has been studied using Monte Carlo simulations and probabilistic analysis. The results of both methods are in a good agreement. Then simple Monte Carlo procedure was employed to study the ordering in products of a polymeranalogous reaction with accelerating effect of neighboring groups. When the reaction with intra- and interchain acceleration and local ordering proceed simultaneously in confined conditions, the ordering might affect the process so that the formation of certain nano-structures (in particular, not trivial strip-like ones) is possible.     

Adsorption of carbon dioxide of 1-site and 3-site models in pillared clays: a Gibbs ensemble Monte Carlo simulation

Adsorption behavior of carbon dioxide confined in pillared clays is analyzed by using constant pressure Gibbs ensemble Monte Carlo (GEMC) method. In our simulation, 1-site and 3-site models are used to represent carbon dioxide. At the 1-site model, carbon dioxide is described as a Lennard-Jones (LJ) sphere, while at the 3-site model, carbon dioxide is modeled as a three-sites linear chain represented by EPM2 potential considering the quadrapolar effect. The potential model from Yi et al. for pillared clays is used to emphasize its quasi two-dimensional structure. Comparing the calculated results from the 1-site and the 3-site models at T=228.15 and 258.15 K, we observe that the adsorption amount from the two models is the same basically. However, the local density presents a significant difference, because the shoulder in the main peak near the wall from 3-site model can reflect the orientation of carbon dioxide. Accordingly, in the systematical investigation to explore the effect of porosity and pore width on the adsorption of carbon dioxide in pillared clays, the 3-site model was only used. We observe that for a narrow pore of H=1.02 nm, each isotherm shape displays type I curve, suggesting that it is not inflected by the porosity. However, for the larger pores of H=1.70 and H=2.38 nm, the increase of the porosity alters the shape of adsorption isotherms from a simple linear relation to the first order jump, indicating that the porosity is of very important factor to affect adsorption and phase behavior of fluids confined in pillared clays. The excess adsorptions of carbon dioxide at supercritical temperatures of T=323.15 and 348.15 K are also investigated. We find that the maximum exists for each excess isotherm, and the optimal pressure corresponding to the maximum increases with the pore width. However, the porosity has no significant effect on the optimal pressure.     

Combined Monte Carlo and molecular dynamics simulation of hydrated 18:0 sphingomyelin-cholesterol lipid bilayers

We have carried out atomic level molecular dynamics and Monte Carlo simulations of hydrated 18:0 sphingomyelin (SM)-cholesterol (CHOL) bilayers at temperatures of 20 and 50 degrees C. The simulated systems each contained 266 SM, 122 CHOL, and 11861 water molecules. Each simulation was run for 10 ns under semi-isotropic pressure boundary conditions. The particle-mesh Ewald method was used for long-range electrostatic interactions. Properties of the systems were calculated over the final 3 ns. We compare the properties of 20 and 50 degrees C bilayer systems with each other, with experimental data, and with experimental and simulated properties of pure SM bilayers and dipalmitoyl phospatidyl choline (DPPC)-CHOL bilayers. The simulations reveal an overall similarity of both systems, despite the 30 degrees C temperature difference which brackets the pure SM main phase transition. The area per molecule, lipid chain order parameter profiles, atom distributions, and electron density profiles are all very similar for the two simulated systems. Consistent with simulations from our lab and others, we find strong intramolecular hydrogen bonding in SM molecules between the phosphate ester oxygen and the hydroxyl hydrogen atoms. We also find that cholesterol hydroxyl groups tend to form hydrogen bonds primarily with SM carbonyl, methyl, and amide moieties and to a lesser extent methyl and hydroxyl oxygens.     

Monte Carlo simulation of polarizable systems: Early rejection scheme for improving the performance of adiabatic nuclear and electronic sampling Monte Carlo simulations

An early rejection scheme for trial moves in adiabatic nuclear and electronic sampling Monte Carlo simulation (ANES-MC) of polarizable intermolecular potential models is presented. The proposed algorithm is based on Swendsen–Wang filter functions for prediction of success or failure of trial moves in Monte Carlo simulations. The goal was to reduce the amount of calculations involved in ANES-MC electronic moves, by foreseeing the success of an attempt before making those moves. The new method was employed in Gibbs ensemble Monte Carlo (GEMC) simulations of the polarizable simple point charge-fluctuating charge (SPC-FQ) model of water. The overall improvement in GEMC depends on the number of swap attempts (transfer molecules between phases) in one Monte Carlo cycle. The proposed method allows this number to increase, enhancing the chemical potential equalization. For a system with 300 SPC-FQ water molecules, for example, the fractions of early rejected transfers were about 0.9998 and 0.9994 at 373 and 423 K, respectively. This means that the transfer moves consume only a very small part of the overall computing effort, making GEMC almost equivalent to a simulation in the canonical ensemble.     

DNA closed nanostructures: a structural and Monte Carlo simulation study

DNA nanoconstructs are obtained in solution by using six unique 42-mer DNA oligonucleotides, whose sequences have been designed to form a pseudohexagonal structure. The required flexibility is provided by the insertion of two non-base-paired thymines in the middle of each sequence that work as flexible hinges and constitute the corners of the nanostructure when formed. We show that hexagonally shaped nanostructures of about 7 nm diameter and their corresponding linear open constructs are formed by self-assembly of the specifically designed linear oligonucleotides. The structural and dynamical characterization of the nanostructure is obtained in situ for the first time by using dynamic light scattering (DLS), a noninvasive method that provides a fast dynamic and structural analysis and allows the characterization of the different synthetic DNA nanoconstructs in solution. A validation of the LS results is obtained through Monte Carlo (MC) simulations and atomic force microscopy (AFM). In particular, a mesoscale molecular model for DNA, developed by Knotts et al., is exploited to perform MC simulations and to obtain information about the conformations as well as the conformational flexibilities of these nanostructures, while AFM provides a very detailed particle analysis that yields an estimation of the particle size and size distribution. The structural features obtained by MC and AFM are in good agreement with DLS, showing that DLS is a fast and reliable tool for characterization of DNA nanostructures in solution.     

Solvation of sodium-chloride ion pair in water cluster at atmospheric conditions: grand canonical ensemble Monte Carlo simulation

Open statistical ensemble simulations are used to study the mechanism of nucleation of atmospheric water on sodium-chloride ion pair in a wide range of temperature and relative humidity values. The extended simple point-charge model is used for water molecules. Ions-water nonadditive interactions are taken into account by introducing the mutual polarization of ions and water in the field of each other. Gibbs free-energy variations are calculated from Na+-Cl- pair-correlation function and used as a criterion for determining the possible stable states of the cluster. In this relation, it was found that the dissociation of ion pairs in water clusters occurs even at vapor pressures of only a few millibars. In the conditions under consideration solvent-separated ion-pair states are found to be more probable than contact ion-pair configurations. The susceptibilities of water and ions are found to play an essential role in the stabilization of ions at large separations. The structure of ion-induced clusters is analyzed in terms of binary correlation functions. The non-pair interactions influence essentially the structure of ion solvation shells. The results of simulation show that the separation of the charges in water clusters containing simple ions can take place under atmospheric conditions.