Parametrization of coarse grained force fields for dynamic property of ethylene glycol oligomers/water binary mixtures
aa r X i v : . [ m a t h - ph ] A ug Parametrization of coarse grained force fields fordynamic property of ethylene glycol oligomers/waterbinary mixtures
Tamio Yamazaki ∗ Analysis technology center, Canon Inc. 3-30-20, Shimomaruko, Ota-ku, Tokyo 146-8501,Japan.
Abstract
To evaluate shear viscosity of ethylene glycol oligomers (EGO)/water binarymixture by means of coarse-grained molecular dynamics (CG-MD) simu-lations, we proposed the self-diffusion-coefficient-based parameterization ofnon-bonded interactions among CG particles. Our parameterization proce-dure consists of three steps: 1) determination of bonded potentials, 2) scal-ing for time and solvent diffusivity , and 3) optimization of Lennard-Jonesparameters to reproduce experimental self-diffusion coefficient and densitydata. With the determined parameters and the scaling relations, we evalu-ated shear viscosities of aqueous solutions of EGOs with various molecularweights and concentrations. Our simulation result are in close agreementwith the experimental data. The largest simulation in this article corre-sponds to a 1.2 µ s atomistic simulation for 100,000 atoms. Our CG modelwith the parameterization scheme for CG particles may be useful to studythe dynamic properties of a liquid which contains relatively low molecularweight polymers or oligomers. Keywords:
Coarse-grained molecular dynamics simulation, Ethylene glycololigomer, Self-diffusion coefficient, Shear viscosity ∗ Corresponding author. Tel.: +81-3-3758-2111.
Email address: [email protected] (Tamio Yamazaki) . Introduction
Aqueous polymer solutions are widely used in industrial and householdapplications. For example, in ink-jet printing, polymer additives are usedto control fixity and osmosis of the ink droplets to the paper. In order toestimate such properties, it is very important to understand the shear vis-cosity and the diffusivity of the polymer solutions. The molecular dynamics(MD) method at atomic level is the most common technique to estimate theshear viscosity of a solution or the self-diffusion coefficient of a component[1, 2]. However, even if a small oligomer of which the number of monomersis around 10, the longest relaxation time (for usual polymer solutions, it willbe the rotational relaxation time of polymer chains), is several tens nanosec-onds, which is about 1000 times longer than the characteristic time of watermolecules. The rotational relaxation time of the polymer chains increases inproportion to some power of its molecular weight.To obtain the ensemble averaged quantities, the sampling time (the time-average period) should be longer than the longest relaxation time of thesystem. The extreme long relaxation time complicates the estimations of theshear viscosity and the self-diffusion coefficient by using all-atom moleculardynamics (AA-MD) of which typical time-steps is femto-second order.In order to extend the accessible time scale, the coarse-grained (CG) tech-niques have been developed for several decades. The CG techniques, inwhich some atomic groups are represented by a single particle, are widelyused for simulations of meso-scale phenomena in soft materials (lipids, sur-factants and polymers), for example, vesicle formation and fusion [3, 4], self-assembly [5, 6], micelle formation [7, 8, 9, 10] and pore formation in lipidbilayer [11, 12, 13, 14]. The coarse grained molecular dynamics (CG-MD)can achieve speed-up of up to several orders of magnitude faster than anAA-MD. Of course, the degree of speed-up depends on the details of the CGmodel or the system size. CG-MD method is very powerful tool to studyof the static property, such as equilibrium structure of a large system, andit is also successful in investigating the fundamental mechanism of the longtime-scale behaviors of soft materials.To apply this method to a quantitative estimation of the dynamics, itshould be noted that the time in CG-MD trajectory is not equivalent to thetime in the all-atom (AA) description, due to the lack of atomistic details inthe CG model. An effective method to obtain the consistency of the molec-ular motions between AA-MD and CG-MD is to multiply the time scale of2he CG-MD by a constant factor.Recently, several authors reported the systematic studies of the dynamicaland rheological properties of polymer systems polymer systems(polystyrenemelt [15, 16], polyamide-6,6 melt [16] and aqueous polyethylene glycol so-lutions [17]) by the multiscale approach that combines AA-MD and CG-MD simulations through the scaling of time in the CG model. With thewell-tuned potential energy function among the particles in the CG model,CG-MD gives us the reasonable possibility to investigate both statical anddynamical properties of large scale systems which include macromolecules.However, there are few studies about the parametrization for the non-bondedCG interactions, which can be adapted to the estimations of dynamic prop-erties of system, especially of the multi-component system.In this article, we will present the results of the CG-MD simulations, in-cluding the systematic determination of the non-bonded parameters for theethylene glycol oligomer (EGO)/water binary mixtures based on the self-diffusion coefficient and density data. We will also present the results ofthe shear viscosity of the mixture estimated by means of the nonequilibriumCG-MD simulations with the newly determined force field parameters.The article is organized as follows. The experimental methods for measure-ment of the shear viscosity and self-diffusion coefficient are explained in thesection two. The explanation of the CG model and the force-field parame-terization for EGO/water binary mixtures and the computational details arein the 3rd section. The results of the comparison between our simulationsand the experimental measurements are shown in the Results and discussionsection that follows.
2. Experimental section
Reagent grade diethylene glycol (DEG), tetraethylene glycol (TEG) andPEG600 were purchased from Sigma Aldrich Chemical Company and usedin this work without further purification. EGO/water binary mixtures wereprepared gravimetrically using distilled water. D O with a minimum deuter-ation degree of 99.95 % (Merck Co. & Inc., Darmstadt, Germany) was usedfor all experiments as the NMR lock solvent.3 .1.2. Experimental details
The shear viscosities of the EGO aqueous solutions were measured at 293K by the RE80 viscometer manufactured by Toki Sangyo Co., Ltd.The self-diffusion coefficient measurements were performed at 293 K by pulsefield gradient spin echo (PGSE)[18] using standard ledbpgp s sequence onthe Avance600 NMR spectrometer (Bruker BioSpin GmbH, Rheinstetten,Germany). To avoid mixing H O in the aqueous EGO solution and D Oas the NMR lock solvent, we used NMR tubes which have a double tubestructure. In the NMR tube, the inner tube was filled with the NMR locksolvent D O, and the outer tube was filled with aqueous EGO solution.In H-NMR spectra of aqueous EGO solutions, the peaks of around δ ≈ . CH O groups of EGO, and the peaksof around δ ≈ . D EGO ( δ ≈ . D OH ( δ ≈ . D OH .Under the assumption that the proton exchange between hydroxyl of EGOend and water molecules is very rapid, a self-diffusion coefficient of water( D w ) can be obtained by [19] D w = χ i − χ i D OH − − χ i D EGO , (1)where χ i is the molar fraction of hydroxylic protons of EGO ends, and D w is the self-diffusion coefficient of water molecule.
3. Theoretical section
In our coarse-grained model, EGO molecules and water molecules arerepresented by coarse-grained particles, as shown in Figure 1.EGO molecules are modeled by two types of particles (”PA”) and (”PB”).The PA particle represents the oxyethylene unit of both ends of a ethyleneglycol chain, and the PB particle represents the oxyethylene unit in a ethy-lene glycol chain. The mass of PA and PB particles are 53 amu and 44 amurespectively. A mass of an oxygen atom of ether is distributed to two ad-joined coarse grained particles at an equal ratio.4 igure 1: Coarse-graining scheme of a tetraethylene glycol and water. Water is modeledas a PW particle which corresponds to three water molecules.
Water is modeled by a single particle (”PW”) corresponding to three real wa-ter molecules. The mass of PW particle is 54 amu. DEG, TEG and PEG600molecules are modeled as PA-PA, PA-(PB) -PA and PA-(PB) -PA, and weexpressed these EGO species as EGO2, EGO4 and EGO13 respectively. We assume that the total potential energy for CG molecule is written as U tot = X U b ( L ij ) + X U ang (Θ ijk ) + X U non ( R ij ) , (2)where the terms, U b , U ang , U non are the effective potential functions of bondlength L ij , bond angle Θ ijk , and the distance between non-bonded CG parti-cles R ij respectively. Here, i , j and k are the indices of the CG particles. Thenon-bonded CG interactions U non ( R ij ) are modeled using the 12-6 Lenard-Jones potential, U non ( R ij ) = 4 ǫ ij (cid:18) σ ij R ij (cid:19) − (cid:18) σ ij R ij (cid:19) ! , (3)where σ ij is the (finite) distance at which the inter-particle potential is zero, ǫ ij is the depth of the potential well. We use the Lorentz-Berthelot combi-nation rule for interaction between the different species, ǫ ij = √ ǫ i · ǫ j , σ ij = σ i + σ j . (4)In the case the interaction between PB and PW particles, ǫ ij is calculatedfrom the eq.(5) instead of eq.(4), ǫ ij = γ · √ ǫ i · ǫ j . (5)5 .3. Scaling relations As mentioned in Introduction, whenever the CG model is used, due tothe lack of atomistic details, the speed of the time-evolution of the CG-MDtrajectory is inconsistent with that of the AA-MD trajectory. To estimateof dynamical properties of the aqueous EGO solutions by means of CG-MD,we introduced the time mapping parameter S , which is defined by Kremeret.al.[15] as a ratio between the effective segment frictions in the AA modeland in the CG model, S ≡ ζ AA ζ CG , (6)where ζ AA is an effective scalar friction coefficient of a segment in the AAmodel, and as ζ CG is ones in the CG model. In this article, we have assumedthat the ratio between the shear viscosities in the AA model η AA and in theCG model η CG is equal to the parameter S , namely: η AA = S · η CG . (7)According to the Stokes-Einstein relation, the self-diffusion coefficient D isinversely proportional to η and the hydrodynamic radius of a segment r h , D ∝ ( η · r h ) − . (8)In our CG model, three water molecules are included in a single PW particle.The PW particle in the CG model has a larger hydrodynamic radius thanthe individual realistic water molecule, and thus the self-diffusion coefficientof water molecule in CG model is somewhat smaller than that in the AAmodel, even if we apply the time mapping parameter S properly.We define the hydrodynamic radius r CGh of the molecule in CG model by r CGh = √ n · r AAh , (9)where n is a number of atomistic molecules, which are included in one coarse-grained molecule, n = 1 (for EGO), n = 3 (for water) as shown in Figure1, and r AAh is the hydrodynamic radius of molecule in the AA model, thefactor √ n in eq.(9) comes from the assumption of the linear relationshipbetween the cube of the hydrodynamic radius of a molecule and its volume.From eqs.(7) and (8), the scaling relation between self-diffusion coefficientsof molecular systems, which are described by the AA model and by the CGmodel, is given by D AA = S − · √ n · D CG , (10)6here D CG is defined using the Einstein relation, D CG = lim t →∞ t D(cid:0) r CGcom ( t ) − r CGcom (0) (cid:1) E , (11)where r CGcom ( t ) is coordinates vector of the center of mass of the molecule inCG model at time t , h· · ·i denotes the ensemble average. In this article, S isdealt with as a constant value for simplification, though the large molecularweight dependence of S is argued in the CG-MD study of polystylene meltin Ref.15.In order to determine S , we consider the self-diffusion coefficient of purewater using CG-MD simulation D CGwater , and the experimental value D EXPwater observed by NMR measurement. From eq.(10), S is given by S = √ · D CGwater D EXPwater . (12) The Boltzmann inversion method is well-known as one of the techniquesto evaluate effective mesoscale potential from atomistic simulation [20]. Thistechnique enables us to determine the bonded intramolecular interactionespecially. In this article, U b ( L ij ) and U ang (Θ ijk ) are determined by thismethod. In a thermal equilibrium system, we assume that an appearanceprobability P of a state vector of a system Q , obeys the Boltzmann distri-bution. P ( Q ) ∝ exp( − βU ( Q )) , (13)where U ( Q ) is a certain effective potential energy as a function of Q , and β is 1 /k B T .Once P ( Q ) is obtained, a effective potential U ( Q ) can be straightforwardlydetermined from the inversion of eq.(13). In this article, we assumed thatall of the bonded (stretching / bending) potentials in the EGO chain areapproximately given by the same U , as we will show in eq.(15).To obtain P ( L ij ) and P (Θ ijk ), the AA-MD simulation of a triethylene glycoldimethyl ether (TEGDE) in gas-phase is performed at 293 K. With ourcoarse-graining manner, a TEGDE molecule is modeled by three coarse-grained particles (PB p , p = 1 , ,
3) as shown in Figure 2.The center of PB particle is defined as the center of mass of the oxyethy-lene monomer unit, in which two oxygen atoms at both ends of the unit are7 igure 2: Coarse-graining scheme of a triethylene glycol dimethyl ether molecule. weighted by 0.5 and other atoms are weighted by 1.0. Histogram H ( L ij ) ofbond length (bond PB − PB and bond PB − PB ) and histogram H (Θ ijk )of bond angle(bond angle PB − PB − PB ) are obtained from the 50-nstrajectory of the atomistic molecular dynamics simulation. Then, these his-tograms are renormalized by P ( L ij ) ∝ H ( L ij ) L ij , P (Θ ijk ) ∝ H (Θ ijk )sin(Θ ijk ) . (14)As a technical subject for the U ( Q ), because there are a lot of noises inthe potential energy function obtained by the Boltzmann inversion scheme,it should be smooth by some appropriate functions. We assumed that theprobability distribution function P ( Q ) can be expressed by the linear com-bination of Gauss functions G l [21]. Then the effective potential energyfunctions U ( L ij ) and U (Θ ijk ) are obtained by U ( Q ) /k B T = − ln( P ( Q )) = − ln( m X l =1 G l ) + const. = − ln m X l =1 A l ξ l p π/ (cid:18) − ( Q − µ l ) ξ l (cid:19) + const. , (15)where A l , µ l and ξ l are total area, center position and width of the Gaussfunction G l respectively, and m is the number of the Gauss functions forsmoothing of P ( Q ). In this article, we decided m = 3 for the both potentialsof bond length ( Q = L ij ) and bond angle ( Q = Θ ijk ). The areas ( A l ), centercoordinates ( µ l ) and widths ( ξ l ) are determined by curve fitting of the cal-culated bond length / bond angle distribution data.8 .5. Parameterization of Non-bonded potentials In our CG model for the EGO / water binary system, there are threedifferent particles (PA, PB, PW), and 7 parameters for the non-bonded in-teraction U non ( R ij ), which are ǫ PW ǫ PA , ǫ PW , σ PW , σ PA , σ PB and γ . Thesenon-bonded parameters for the U non ( R ij ) are determined based on the ex-perimental data (the densities and the self-diffusion coefficients for severalEGO / water binary systems). The parameterization procedure consists ofthe four steps listed below.STEP 1 : The parameter ǫ PW as a function of parameter σ PW is determinedso as to reproduce the density of pure water at 293 K, which is 0.998 g / cm [22]. In order to satisfy the experimental pure water density, ǫ PW is uniquelydetermined according to σ PW . The self-diffusion coefficient of pure water D CGwater is evaluated by CG-MD simulation with the fixed parameter values( ǫ PW and σ PW ). The time mapping parameter S is calculated from eq.(12).STEP 2 : The parameter ǫ PA as a function of parameter σ PA is determinedso as to reproduce the density of pure diethylene glycol (EGO2) liquid at 293K, which is 1.118 g / cm [23]. With the parameters for the PW determinedin STEP 1, unique σ PA is obtained through minimizing of the error betweenthe D CG calculated by CG-MD and the D AA observed by NMR of the com-ponents in the EGO2/water binary mixtures (EGO2 weight fraction : 0.2,0.5 and 0.8).STEP 3 : For determination of the parameter ǫ PB , we used the density ofpure tetraethylene glycol (EGO4) liquid at 293 K, which is 1.125 g / cm [23]and the self-diffusion coefficients of the components in the triethylene glycol(EGO3)/water binary mixtures (EGO3 weight fraction : 0.2, 0.5 and 0.8)observed by NMR. The EGO4 and EGO3 are modeled in this article as PA-PB-PB-PA and PA-PB-PA, respectively. With the parameters for the PWand the PA determined already in STEP 1 and STEP 2, the σ PB can beobtained by the same process in STEP 2. At the first execution of STEP3, γ of eq.(5) is setted to 1. After performing STEP 4, γ is revised to theoptimized value.STEP 4: In the case of the calculating of the interaction between PB andPW particle, the parameter γ , which adjust the miscibility of the EGO chainin water, should be determined. The γ is obtained through minimizing ofthe error between the D CG calculated by CG-MD and the D AA observed byNMR of the components in the PEG600 (EGO13)/water binary mixtures(EGO13 weight fraction : 0.2). 9TEP 3 and 4 are sequentially repeated until the proper σ PB and γ areobtained. In STEP 2, 3 and 4, the root means of square errors (RMSE) areevaluated by RM SE = vuut M M X q =1 (cid:0) log D EXP q − log D AA q (cid:1) , (16)where M is the number of data, D EXP q is experimental self-diffusion coeffi-cients, D AA q is calculated self-diffusion coefficients from eq.(10). All simulations of this work are performed by GROMACS 4.0.5 [24]. Inthe atomistic molecular dynamics simulation of gas phase of TEGDE, thetemperature is held at T = 293 K by Langevin thermostat [25] with cor-relation time τ = 1.0 ps . The production time step for integration is dt= 1 fs. And cutoff radius for LJ and Coulomb potentials is 1.4 nm. Thegeneral Amber force field (GAFF) [26] is used as the atomistic force fieldsfor TEGDE molecule. Atomic charges of a TEGDE molecule are assignedby AM1-BCC [27] method. These force fields and atomic charges are gen-erated by antechamber 1.4 [28, 26] and acpype 1.0 [29]. In the CG-MDsimulations, to evaluate the density and the equilibrated liquid structure, 2ns MD simulations are performed in the N pT ensemble. Nose-Hoover ther-mostat [30, 31, 32] is used at 293 K to control temperature of the system.Parrinello-Rahman barostat [33] is used at 1 atm to control pressure of thesystem. In the 2 ns MD simulation, the instantaneous densities are calcu-lated from the last 1 ns trajectory every 1000 steps and then are averaged.Omitting the first 1 ns data of trajectory as relaxation time, the self-diffusioncoefficient is calculated from 3 ns (for water, EGO2/water, EGO4/water bi-nary mixture) or 30 ns (for EGO13/water binary mixture) MD simulation,which is performed at constant volume
N V T ensemble, at the initial struc-ture is the last record of the former 2 ns MD simulation.Non-equilibrium molecular dynamics(NEMD) [34] simulation is applied forthe calculation of the shear viscosity of aqueous EGO solution. In order tocalculate the shear viscosity, 20 ns (for water, EGO2/water, EGO4/waterbinary mixture) or 200 ns (for EGO13/water binary mixture) MD simula-tion is performed. After dropped first 5 ns trajectory, the shear viscosity iscalculated by analysis of NEMD trajectory.10here are 5 sets of MD/NEMD simulations with different initial structuresand initial velocity profiles. The diffusion coefficients or the shear viscosi-ties are calculated from each MD/NEMD simulations, and then are aver-aged. The cutoff radius R cut = 1.4 nm for non-bonded interaction amongthe coarse-grained particles and a production time step dt = 10 fs for inte-gration of Newton’s equation are used as common conditions in all CG-MDsimulations.
4. Results and discussion
Parameters ( A l , µ l and ξ l ) of eq.(15), which are determined from curvefitting of the data of U ( L ij ) and U (Θ ijk ) are summarized in Table 1 andTable 2, respectively. Table 1: Parameters of bond length potential represented by eq.(15) bond type l A l µ l (nm) ξ l (nm)PA-PA,PA-PB and PB-PB 1 0.382 0.023 0.3052 0.229 0.020 0.3383 0.028 0.018 0.266 Table 2: Parameters of bond angle potential represented by eq.(15) angle type l A l µ l ( ◦ ) ξ l ( ◦ )PA-PB-PB and PB-PB-PB 1 238.840 57.471 190.5672 45.375 24.819 123.9863 31.826 14.765 101.560 Figure 3 shows the comparison between the AA model and the CG modelfor the histograms of the length L ij (a) and of the angle Θ ijk (b) of bondedPB particles of TEDME, as shown in Figure 2. The Lennard-Jones parameters for PW, PA and PB particles, the pa-rameter γ of eq.(5) and the time mapping parameter S are straightforwardlydetermined by using the procedures (STEP 1 to STEP 4) mentioned as sec-tion 3.5. In STEP 1, to satisfy the experimental liquid density of pure water11 igure 3: Probability distributions of bond length between bonded PB particles (a) and ofbond angle among bonded PB -PB -PB particles (b) of TEDME. Open Symbols denotethe distribution obtained by using AA-MD simulation. Broken line denotes the distribu-tion obtained by CG-MD simulation. (0.998 g/cm ) at 293 K and 1 atm pressure, ǫ PW is uniquely determinedaccording to σ PW , which is limited in range from 0.38 to 0.42 nm. If theparameter σ PW is within this range, then the PW fluid has a stable liquidphase. Though the σ PW can be selected arbitrary within this range, we se-lected σ PW = 0.40 nm. The observed self-diffusion coefficient of pure water at273 K and 1 atm pressure is D EXPwater = 2.0 m /s by using NMR. From eq.(12), S = 6.19 is obtained from eq.(12). Through STEP 2 to STEP 4, the σ PA , σ PB and γ parameters are one by one determined by minimizing the RMSE,which evaluated by eq.(16). These optimized non-bonded parameters aresummarized in TABLE 3. 12 able 3: Parameters of non-bonded potential represented by eqs.(3)-(5) σ P W σ P A σ P B ǫ P W ǫ P A ǫ P B γ From eq.(11), the self-diffusion coefficient is obtained by evaluating ofthe slope of the mean square displacement (MSD) of the center of mass ofthe EGO/water molecules. The MSDs were calculated for each EGO andwater molecules in the simulation box using 3 ns (for water, aqueous EGO2and EGO4 solutions) or 30 ns (for aqueous EGO13 solution) MD-simulation.Dropping the first 1 ns MSD data as relaxation time, the slope of MSD versustime was evaluated by linear regression.Figure 4a and 4b show the comparison between experimental and calculatedself-diffusion coefficients. Figure 4a shows the self-diffusion coefficients ofthe EGO molecules plotted against the EGO weight fraction ( W ) of theEGO/water binary mixtures. Figure 4b shows the self-diffusion coefficientsof the water molecules plotted in the same manner as Figure 4a. Three dif-ferent EGOs, which are EGO2, EGO4 and EGO13, and three different EGOweight fractions ( W = 0.2, 0.5 and 0.8) are presented, including the data ofaqueous EGO13 solutions ( W = 0.5 and 0.8), which did not used in theparametrization of non-bonded potentials at Section 4.2.Experimental observations show that the self-diffusion coefficients (both ofthe EGO and the water) are linear on a log scale when plotted against W .The CG-MD simulation results are also linear as shown in Figure 4a and4b. The dependence of the D on W is increasing against the EGO molec-ular weight ( M W ) as shown in Figure 4a, but the self-diffusion coefficientsof the water molecule does not show such tendency (Figure 4b). The ex-perimental W and M W -dependences of D are reproduced by the CG-MDsimulations correctly. Additionally, we performed CG-MD simulation ofEGO45(PEG2000)/water binary mixture ( W = 0.2). We found that thecalculated self-diffusion coefficients of EGO45 and water are comparable tothe experimental values based on NMR [19] (data not shown), though theseNMR measurements were observed at 298 K, differ from our CG-MD condi-tion T = 293 K.The shear viscosity was calculated by the nonequilibrium method described13 igure 4: Self-diffusion coefficients D of EGO (a) and water (b) molecules as a functionof the weight fraction W of EGO in the EGO/water binary mixtures. The open circles,the open triangles and the open squares denote the D EXP values measured using NMRfor EGO2/water, EGO4/water and EGO13/water binary mixtures, respectively. Thefilled circles, the filled triangles and filled squares denote the D CG values calculated usingCG-MD simulations for EGO2/water, EGO4/water and EGO13/water binary mixtures,respectively. by B. Hess [34], which estimates the shear viscosity of liquid from nonequi-librium simulation with an external cosine acceleration profile of the form a x (z) = A cos ( k z) , (17) k = 2 πl z , where a x is the acceleration in the x direction, l z is the height of the simulationbox, A is the magnitude of acceleration, z is the coordinate z-direction. The14hear viscosity of liquid η can be estimated by η − ( t )= 2 k ρA N X i =1 m i v i, x ( t )cos( k z i ( t )) / N X i =1 m i , (18)where ρ is the density of a system, m i is the mass of the particle of index i , v i, x is the velocity in the x direction of the particle of index i , N is the totalnumber of particles in the simulation box, η − is the reciprocal of viscosityof the system. The A parameter must be carefully selected, the shear rateshould not be so high that the system is driven too far from equilibrium.The maximum share rate of the CG system sh CGmax issh
CGmax = A ρη CG l z π . (19)For our simulations with: η CG ≈ kgm − s − ], l z ≈ ρ ≈ − ], and sh CGmax is approximately 0.2 [ps nm − ] A . This shear rate shouldbe smaller than one over the longest correlation time in the system. For usu-ally liquids, it will be the rotational correlation time of the largest molecule inthe system. In the aqueous EGO13 solution ( W = 0.8), the rotational relax-ation time of end-to-end vector of the coarse-grained EGO13 is approximately6000 ps. In this case, parameter A should be smaller than 0.001 [nm ps − ].When the shear rate is too high, the observed shear viscosity will be too low.In this article, we used A = 0.0005 [nm ps − ] for all nonequilibrium CG-MDsimulations. η AA is estimated by eq.(9), with S = 6.19 and η CG obtainedfrom eq.(18). Figure 5 shows the comparison between experimental and cal-culated shear viscosities of the aqueous EGO solutions. For nine EGO/waterbinary mixtures evaluated in this article, the symbols are same as in Fig.4,the calculated shear viscosities are agreed with the experiments. This meansthat the calculated values for two EGO-concentrations( W = 0.2 and 0.5)are very close to the experimental data and overall data included W = 0.8shows the same tendency as the experimental results. In order to verify ourCG model at larger molecular weight of EGO, we performed the additionalcalculations of the EGO22(PEG1000)/water and EGO45(PEG2000)/waterbinary mixtures. At the highest EGO-concentrations ( W = 0.5 for EGO22and W = 0.4 for EGO45) considered in this article, the rotational relaxationtimes of end-to-end vector of both the coarse-grained EGO22 and EGO45are lower than the ones of the coarse-grained EGO13 in the aqueous EGO1315olution ( W = 0.8). Therefore, the parameter A in the eq.(18) and the sam-pling time length of η − are same as the case of the EGO13/water binarymixture ( W = 0.8 ). Figure 6 shows the comparison between our calculatedresults and the experimental observations measured at 293 K [35]. We foundthat the calculated viscosities are comparable to the literature values foraqueous EGO22 solutions. Although our calculated viscosities are slightlylower than the literature values for aqueous EGO45 solutions, the tendencyof the dependence of viscosity on the EGO45 concentration seems to be inaccordance with the experiments. Figure 5: Shear viscosity η of the EGO and water binary mixture as a function of theEGO weight fraction W . The symbols are same as in Fig.4. igure 6: Shear viscosity η of the EGO and water binary mixture as a function of the EGOweight fraction W . The open circles and the open triangles denote the η EXP values in lit-erature [35] for EGO22(PEG1000)/water and EGO45(PEG2000)/water binary mixtures,respectively. The filled circles and filled triangles denote the η CG values calculated usingCG-MD simulations for EGO22/water and EGO45/water binary mixtures, respectively.
5. Conclusions
CG-MD simulations for the EGO/water binary mixtures were performedat 293 K and 1 atm pressure. The EGO chain was modeled by two typesof particles, PA and PB, which represent the oxyethylene units of both endsand the middle of the chain, respectively. Also, three water molecules areincluded into a single PW particle. With our CG model, the number of par-ticles that should be considered in the simulation can be reduced ten-fold,and the time step for integration of Newton’s equation increases ten timescompared with atomistic simulations.The parameters for the CG model were determined by the systematic man-ner, as was shown in this article. The CG bonded potential parameters forthe EGO chain were obtained by the Boltzmann inversion of the correspond-ing atomistic distribution functions. Due to the soft pair potential amongthe CG particles, the time-scale of CG-MD simulation is not equivalent to arealistic time scale. In order to estimate the proper dynamical properties bymeans of CG-MD simulations, the scaling relation for time in the CG modelis introduced. Moreover, the hydrodynamic radius of PW particles in ourCG model is larger than those of atomistic water molecules, due to the gainin volume of the PW particle as the result of the coarse graining of water,17s shown in Figure 1. Therefore, we also introduced the scaling relation forthe water diffusivity based on the Stokes-Einstein law.With the bonded potential parameters and the scaling relations, the 12-6Lennard-Jones non-bonded potential parameters are straight-forwardly de-termined to reproduce the experimental observations (the density and theself-diffusion coefficient). We adopted the Lorentz-Berthelot mixing rule forthe non-bonded interactions between the unlike particles. In this article, theLorentz-Berthelot mixing rule is slightly modified, as shown in eq.(4) andeq.(5), to estimate the proper miscibility of the EGO13 in water. By usingthe determined CG force-field parameters for the EGO/water binary mix-tures, our CG-MD simulation gives the estimations which agreed well withthe experimental shear-viscosity data, including of the PEG1000/water andthe PEG2000/water binary mixtures which were not used in our parame-terization procedure. The largest simulation in this article corresponds toa 1.2 µ s atomistic simulation for 100,000 atoms. Our CG model with theparameterization scheme for the CG particles may be useful to study of thedynamic properties of a liquid which contains relatively low molecular weightpolymers or oligomers.In our future work, we plan to investigate the CG models for the sev-eral watersoluble polymers (e.g., polypropylene glycol, polyvinyl pyrrolidone,polyvinyl alcohol), for the estimations of the shear-viscosity and the diffusiv-ity of these water mixtures through CG-MD simulations.18 eferenceseferences