Directly Modifying the Nonbonded Potential Based on the Standard

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Directly Modifying the Nonbonded Potential Based upon the Standard Iterative-Boltzmann-Inversion Method for Coarse-Grained Force Fields Zhimin Xie, Dongliang Chai, Youshan Wang, and Huifeng Tan J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b06457 • Publication Date (Web): 21 Oct 2016 Downloaded from http://pubs.acs.org on October 22, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Directly Modifying the Nonbonded Potential Based upon the Standard Iterative-Boltzmann-Inversion Method for Coarse-Grained Force Fields Zhimin Xie,*,† Dongliang Chai,† Youshan Wang, and Huifeng Tan National Key Laboratory of Science and Technology on Advanced Composites in Special Environment, Center for Composite Materials and Structures, Harbin Institute of Technology, Harbin 150001, People’s Republic of China.

ABSTRACT: The effective potentials are of great importance for the coarse-grained (CG) simulations, which can be obtained by the structure-based iterative Boltzmann inversion (IBI) method. However, the standard IBI method is incapable of keeping the mechanical and thermodynamic properties of the CG model in consistence with those of the all-atom model. Unlike the existed techniques, such as introducing the friction force as the dissipative force to drop the superatom motion while the conservative force arising from the CG potential was kept intact, we directly modified the standard IBI nonbonded potential by adding an empirical function. According to an analysis of the dissipative particle dynamics, the additional function did compensate for the friction reduction of the standard IBI CG model. In the present work, the thermal fluctuation information from the nonbonded radial distribution function was incorporated into the additional empirical function. As an illustration of the new CG force fields,

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 39

we presented simulations of stress-strain relation and thermodynamic properties in terms of cispolyisoprene, and compared the statistical structure information of the superatoms with those of the IBI CG model and the all-atom model. It should be emphasized that the additional empirical function contributed to compensation for the friction reduction no matter what functional form it took. In this sense, the proposed method was easily operational.

1. INTRODUCTION The all-atom simulations remain expensive burdens to calculate the macromolecular properties, because the systematic simulations of huge number of atoms have to be performed with a time step small enough to accurately record the high-frequency vibration of atoms in the condensed phase.1-2 However, the structure and dynamics of many polymers, and biological processes often span a large range of length and time scales.3-7 For the long spatiotemporal phenomena of soft matters, the all-atom computational simulations with high resolution are not always economic to discover the underlying mechanism of associated issues. Useful techniques of surmounting this problem involve the coarse-grained (CG) methods, which are very popular in the simulations of soft matter and biological systems.8-13 The main philosophy behind the most coarse-graining methods is to reduce the structural and interaction information of original atomistic model by clustering a group of atoms into new CG sites, or namely superatoms, with lower resolution. The accessible spatiotemporal scales of CG models bridge the gap between the microscopic and mesoscopic simulations.14


2

Page 3 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The effective potentials are of great importance for the CG simulations, so there has been considerable interest in the study of the CG potentials in recent years.11, 15 The structure-based iterative Boltzmann inversion (IBI)16-17 method is one of the most important ways to derive an effective pair potential, for it provides a robust engine that is automatic and iterative. And the standard IBI method has made success in simulating the structural and some thermodynamic properties of many typical polymers.11, 18-20 In practice, the IBI CG potential is a potential of mean force or a free energy. Rigorously, using the free energy as the potential energy in the IBI CG simulation is incorrect in a viewpoint of the statistical mechanics. The IBI CG potential was generally separated into the bonded potential and nonbonded potential in the simulations.11 The nonbonded interaction potential can be obtained by Boltzmann-inverting the radial distribution function (RDF). As stated by the Henderson’s theorem,21 the difference between two different pair potentials that produce the identical RDFs can be expressed only by a constant; nevertheless, two remarkably different nonbonded pair potentials derived from the standard IBI method can yield almost the same two RDFs.22-23 For the physical problems, the standard IBI method inserts a flexible linear tail function to the nonbonded potential so that the target pressure is reproduced under the condition of little altering the target RDFs.17 When a group of atoms is lumped into and considered as a quasi-particle superatom, the excluded volume and the surface roughness of each molecule are changed,24-25 this more or less results in the abnormal dynamic behavior of the CG systems. In the practical simulations of using the standard IBI potentials, the CG model demonstrated much faster movement of the superatoms1, 25 and much lower elastic modulus26-27 than the parent atomistic model. This meant that the CG models reduced the internal friction force.24 Therefore, to increase the friction between the superatoms may improve the IBI CG simulations. Up to now, many efforts have


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 39

been devoted to the problem of the abnormal dynamic and mechanical behavior of the IBI CG systems, and the effective strategies mainly included the method of time mapping,2 the hybrid method of combining the bottom–up and top–down techniques,9 and the method of introducing the dissipative particle dynamics (DPD).24, 28 The time mapping method assumed that all critical dynamic processes were speeded up by the identical factor.1-2, 29 Then, a simply rescaling of time by a proper time-scale factor can reproduce the same dynamics as the parent atomistic model.1 The single time-scale factor was generally obtained from comparing the dynamical characteristic parameter, such as the friction coefficient or the mean-squared displacement, of the CG model with that of the atomistic model.29-30 The time mapping method has been widely applied to reproduce reasonable dynamic motion for low degrees of coarse-graining.2, hybrid method of combining the bottom–up and top–down techniques,9,

11, 18, 29-38

39-40

In the

some structural

accuracy was sacrificed for the CG model to reproduce the target dynamic or mechanical properties of interest. The typical models were designed by Keten’s research group,39 in which, keeping the IBI bonded potentials (that belong to the bottom-up technique) unchanged, the IBI nonbonded potential was replaced by the Lennard-Jones potential (that belongs to the top–down technique) with two unknown coefficients. The hybrid method was robust and effective to predict the thermal and mechanical properties of such polymers as methacrylate-based polymers and polystyrene.9, 39-40 Additionally, the Kirkwood-Buff IBI method proposed by Ganguly et al.22, 41

and the multistate IBI method as reported by Moore et al.23, 42 were also modifications of the

nonbonded potentials to increase the accuracy of the standard IBI method, through increasing thermodynamic consistency22,

41

and through reducing structural artifacts and allowing

density/pressure relationships to be implicitly captured,23, 42 respectively; however, neither of these modifications have been shown to enhance mechanical properties. It is also useful for the


4

Page 5 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


standard IBI CG simulations to regulate the molecular motion by introducing the DPD thermostat.24 In this technique, the conservative force arising from the IBI CG potential was kept intact, and the friction force was introduced as the dissipative force to drop the velocity of the superatoms, accompanied by a random force as the noise driving force. According to the fluctuation-dissipation theorem, the dissipative force and the random force were interrelated together to serve as a thermostat.43 In the standard IBI algorithm, the resultant force acting on individual particles is the conservative force, while in the DPD algorithm the resultant force includes the dissipative and the random forces besides the conservative force. The difference between the two algorithms results in the difference of the nonbonded RDF; in other words, the nonbonded RDF from the DPD algorithm can be considered as the combination of the nonbonded RDF from the standard IBI CG model with a perturbation. Motivated by the modifications of the nonbonded potentials22, 39, 42

and the DPD algorithm,24, 44 the effect of the dissipative-random forces might be introduced

into the nonbonded RDF in a perturbation form to replicate the effect of the DPD thermostat so as to gain a mechanical consistency without utilizing the DPD thermostat. In this paper, we directly modified the standard IBI nonbonded potential by adding an empirical function as a perturbation to the nonbonded RDF to compensate for the friction reduction of the standard IBI CG model. In this way, the CG force field was easily obtained. As an example, the thermal fluctuation information from the nonbonded radial distribution function was incorporated into the additional empirical function. The molecular dynamics simulations in terms of poly(cis-1,4isoprene) (PI) were carried out to demonstrate the proposed method. 2. METHODS AND EXPERIMENTAL SECTION


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 39

2.1. Standard IBI Method. It is necessary to review the standard IBI method17 briefly. The standard IBI method17 is an important approach of deriving the force fields of the CG models. Figure 1 shows a schematic procedure of the standard IBI method to develop the bonded and nonbonded potentials. Initial setup of coarsegrained(CG) potentials CG simulation

CG probability distributions

Update CG potentials of superatoms No

In consistence with target probability distributions Yes Finial CG potentials Figure 1. Schematic procedure of the standard iterative Boltzmann inversion method. In the IBI method, the CG potentials are expressed as the functions of the structural probability distributions, which include the bond length probability Pbond (l ) , the bond angle probability

Pangle (θ ) , and bond dihedral probability Pdihedral (φ ) . The target structural distributions are commonly used to set up the initial CG potentials which are given by 0 U bond (l ) = − k BT ln( Pbond,target (l ) / l 2 ) ,

(1)

0 U angle (θ ) = −kBT ln( Pangle,target (θ ) / sin θ ) , (2) 0 U dihedral (φ ) = −kBT ln Pdihedral,target (φ ) ,

(3)

0 U nonbonded ( r ) = − k BT ln g target ( r ) ,

(4)


6

Page 7 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in which k B and T are the Boltzmann constant and the absolute temperature, respectively; the subscript ‘target’ refers to from the all-atom simulation; bond-length stretching potential Ubond, bond-angle bending potential Uangle, bond-dihedral torsion potential Udihedral are the bonded potentials, and Unonbonded stands for the CG nonbonded potential. In the beginning, the initial probabilities of the CG simulation are assumed to be the target ones. From the initial CG simulation, however, the first statistical results of the probability distributions often differ from the target ones. Therefore, the improvement for CG potentials is needed, and this can be updated by adding correction terms from the probability comparison of the current CG simulation with the target model. The improved potentials are used in the next CG simulation step. Similarly, this process of improving the CG potentials is repeated and iterated as follows:

U

n +1 bond

(l ) = U

n bond

n Pbond (l ) / l 2 (l )+kBT ln ， (5) Pbond,target (l ) n Pangle (θ )

n +1 n U angle (θ ) = U angle (θ )+kBT ln

(6) sin θ ， Pangle,target (θ )

n +1 n U dihedral (φ ) = U dihedral (φ )+kBT ln

U

n +1 nonbonded

(r ) = U

n nonbonded

n Pdihedral (φ ) ，(7) Pdihedral,target (φ )

g n (r ) (r )+kBT ln ，(8) g target (r )

in which the superscript ‘n’ denotes the iteration step. The nonbonded potential of the superatoms can be further rewritten as, n g (r )   U nonbonded (r ) = − k BT ln  g target (r ) ⋅ ∏ target  . (9) i g (r )  i =1 


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 39

The iteration process continues until the target distributions from the all-atom model are reproduced by the CG simulation. And the final standard IBI CG potentials were then obtained after the pressure correlation.17, 45

2.2. Modification of the Nonbonded Potential Function. The standard IBI potentials were primarily utilized to compute the conservative force, but ignored the contributions of the dissipative force and random force between the particles. The velocity-Verlet algorithm46 for the standard IBI CG model is,

f (t ) 1 ri (t + ∆t ) = ri (t ) + ∆tvi (t ) + (∆t )2 i ，(10) 2 mi

1 f (t ) + fi (t + ∆t ) v i (t + ∆t ) = v i (t ) + ∆t i ，(11) mi 2 where r and v are position and velocity vectors, respectively; t is the current time; ∆t is a small

(

)

time step; m stands for the mass of superatom; fi (t ) = ∑ FijC (t ) is the resultant force acting on j ≠i

individual particles, and FijC (t ) is the conservative force between particle i and particle j within the cutoff radius rcutoff. The RDFs of the i-th superatom at the beginning moment t and at the ending moment t + ∆t are described by g i ( r , t ) and g i ( r , t + ∆t ) , respectively. According to the studies of Qian et al.,24 the incorporation of the contributions of the dissipative and random forces as a thermostat condition into the physical property calculations can effectively slow down the abnormal acceleration of the individual superatoms in the standard IBI CG simulation, and then a satisfied result was obtained. Specifically, the resultant force in the DPD algorithm44, 47 includes the dissipative FijD (t ) and the random forces FijR (t ) besides the


8

Page 9 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(

)

conservative force, i.e., fi (t ) = ∑ FijC (t ) + FijD (t ) + FijR (t ) . The three types of forces within the j ≠i

cutoff radius rcutoff were expressed as follows,44, 48

FijC = Fij( c ) (rij )eij ,

(12)

FijD = −γ wD (rij )(eij ⋅ vij )eij ,

(13)

FijR = σ wR (rij )ξij ∆t −1/2eij ,

(14)

where rij = ri − r j , rij = rij , eij = rij / rij , v ij = v i − v j ; Fij( c ) (rij ) is specified by the standard IBI

CG potential in this work. The variables γ and σ are the strengths of the dissipative and random forces, respectively, and σ 2 = 2γ k BT ; ξij is a Gaussian random number with zero mean and unit variance, which is independent for each pair of particles and at every time step; wD (rij ) and 2

wR (rij ) are weighting functions, satisfying wD (rij ) =  wR (rij )  and wR (rij ) = 1 − rij / rcutoff .44 A modified version of the velocity-Verlet algorithm44 for the DPD is,

riDPD (t + ∆t ) = riDPD (t ) + ∆tv iDPD (t ) 1 + ( ∆t ) 2 fiDPD (t ) / mi , 2

(15)

v% iDPD (t + ∆t ) = viDPD (t ) + λ∆tfiDPD (t ) / mi , (16) fiDPD (t + ∆t ) = fiDPD (r DPD (t + ∆t ), v% DPD (t + ∆t ))， (17)

1 v iDPD (t + ∆t ) = v iDPD (t ) + ∆t fiDPD (t ) 2 DPD +fi (t + ∆t )  / mi ,

(18)

in which fiDPD (t ) is the resultant of the three types of forces of the i-th superatom. An empirical integration constant λ is used to control the velocity change of the i-th superatom.


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 39

Therefore, the effect of the DPD thermostat can be attributed to the difference between the two algorithms. To capture this difference, we set the initial conditions of the two algorithms in full consistence with each other, e.g., riDPD (t ) = ri (t ) , viDPD = vi (t ) , and giDPD (r , t ) = gi (r , t ) . Additionally, both of the conservative forces in the two algorithms were derived from the same standard IBI potentials. At the next time t + ∆t , the position difference of the i-th particle between the DPD and the standard IBI CG algorithms yields, ∆ri (t + ∆t ) = riDPD (t + ∆t ) − ri (t + ∆t ) f DR (t ) 1 = (∆t )2 i , 2 mi

(19)

in which fiDR (t )=∑ ( FijD (t ) + FijR (t ) ) . ∆ri (t + ∆t ) is bound to result in the nonbonded RDF j ≠i

difference at the time t + ∆t between the DPD and the standard IBI algorithms. Thus, the nonbonded RDF g DPD of the DPD system at the time t + ∆t can be described as,

g DPD (r , t + ∆t ) = g (r , t + ∆t ) − ∆g (r , t + ∆t )，(20) Over a long enough period of time, the more general relationship between g DPD ( r , t + ∆t ) and g ( r , t + ∆ t ) for one time step ∆t can be averaged, i.e.,

g DPD (r ) = g (r ) − ∆g (r , ∆t )  ∆g (r , ∆t )  = g (r ) 1 − ， g (r )  

(21)

where g DPD ( r ) and g ( r ) are the time-averaged nonbonded RDFs from the DPD and IBI CG algorithms at any moment, respectively. Similar to eq 9, the nonbonded potential from the DPD algorithm can be described as, n g  target ( r )  DPD U nonbonded (r ) = −k BT ln  g target (r )∏ l ,DPD  . (r )  l =1 g 

(22)


10

Page 11 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∞

Using eq 21 and Taylor series formula 1/ (1 − x ) = ∑ x n， ( x < 1) , eq 22 is further rewritten as, n =0

 DPD U nonbonded (r )= − k BT ln  g target (r )  n g ∆g l (r , ∆t )   target ( r )  ⋅∏ l 1 +  . g (r )  g l (r )   l =1 (23) Let

n

g target (r )

i =1

g i (r )

g IBI ( r ) =g target (r ) ⋅ ∏

, defined as the effective IBI RDF. The term

n

∏ (1 + ∆g (r , ∆t ) / g (r ) ) l

l

in eq 23 can be represented by the equivalent expression

l =1

1 + G wave ( r , ∆t ) / g IBI ( r ) , which means that the additional function G wave ( r , ∆t ) is directly related to the coupled dissipative-random forces just like ∆g l ( r , ∆t ) . After introduction of 1 + G wave ( r , ∆t ) / g IBI ( r ) , the eq 23 is transformed into DPD U nonbonded (r )

 G wave ( r , ∆t )   = − k BT ln  g IBI ( r ) ⋅ 1 +  g IBI ( r )    = − k BT ln  g IBI ( r ) + G wave ( r , ∆t )  .

(24)

DPD As clearly seen through a comparison of eqs 24 and 9, the nonbonded potential U nonbonded (r )

from the DPD algorithm can be reduced to that from the standard IBI algorithm when

G wave ( r , ∆t ) = 0 . So, eq 24 indicates that the additional perturbation G wave ( r, ∆t ) can be directly introduced into the nonbonded force field to compensate for the friction reduction of the standard IBI CG model, no matter what functional form it took. In this sense, the proposed method was easily operational. It is anticipated that the direct modification of the standard IBI nonbonded


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 39

potential would recover more than one essential property of the target all-atom model, such as the mechanical property in addition to the pressure. To some extent, the motion of a superatom under complex nonbonded interactions can be treated as a kind of Brownian motion, which is governed by the fluctuation-dissipation theorem. Thanks to the fluctuation-dissipation theorem,49 the dissipative and random forces as the external perturbation between the local nonbonded superatoms lead the nonbonded superatom into the fluctuating state with thermodynamic equilibrium. Therefore, in the present work, the thermal fluctuation information from the nonbonded RDFs is incorporated into the additional empirical function. As an example, the additional function G wave ( r, ∆t ) was taken as a sinusoidal function, which was abstracted from the primary characteristics of the nonbonded RDF fluctuations induced by the local friction effects. This simplification of the nonbonded RDF fluctuations mainly took three indispensable aspects into account to keep the crucial physical information. The first key aspect was that the thermal fluctuations show some spatiotemporal periodicity, and therefore, a sinusoidal function as the main body of G wave ( r, ∆t ) was a good choice to achieve this periodic functionality. The second aspect was about the fluctuation amplitude of nonbonded RDF, which was important for its dominating the magnitude of damping interactions between the nonbonded superatoms. Based on the fitted result of the amplitude outline, a compound function that included two adjustable parameters was to be constructed to optimize the whole sinusoidal function of G wave ( r, ∆t ) . The third aspect was related to the constrain condition. Since a higher value of g(r) corresponded to more frequent and strong RDF fluctuations for the nonbonded superatoms that had a lower potential, i.e., the period and the amplitude of the thermal fluctuations of the nonbonded RDF were also directly proportional to the value of g(r), the


12

Page 13 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


period and the amplitude of the sinusoidal function should be zero when g(r) is equal to zero at an extreme condition. Integrating all the aforementioned aspects, the additional empirical function of nonbonded RDF in our present work was given by

G wave ( r , ∆t ) ≈ g 0 (r ) f outline ( r )  2π  ⋅ sin  ( r − rfirst ) + ϕ0  , η ( r ) rcycle 

(25)

in which f outline ( r ) was an adjustable compound function based on the amplitude outline controlling the amplitude of the thermal fluctuations of nonbonded RDF; η ( r ) was a probability ratio of g 0 (r ) to g 0 ( rfirst ) ; rfirst was the corresponding distance r of the first peak value of g 0 (r ) ;

rcycle and ϕ 0 were the cycle distance for controlling the period and the phrase angle of the sinusoidal function at the distance of rfirst , respectively. The additional empirical function

G wave ( r, ∆t ) in eq 24 may also take other physical forms, which was not unique. Additionally, the variation of the nonbonded RDF will change the pressure of the CG system; hence the IBI pressure-correction method17,

44

was used to recover the pressure at any time of nonbonded

potential modification. 2.3. Description of Simulation. The all-atom model consisting of 100 chains of PI with 10

monomers in every chain was constructed in the 54 Å × 54 Å × 54 Å periodic boxes by using the Amorphous Cell Module of Materials Studio. The all-atom simulation was carried out by Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) using a consistent valence force field (CVFF)50 with a cutoff radius rcutoff of 15 Å. The temperature-rising/annealing procedure was adopted to accelerate the relaxation process of the system with 1 fs as the time


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 39

step. Following the minimization, the all-atom model was at first equilibrated for 5 ns under NPT ensemble at the higher temperature of 598 K and a pressure of 1 atm; and then, keeping the constant pressure, the all-atom model was relaxed for the next 10 ns under NPT ensemble at the lower temperature of 298 K to achieve a new equilibrium. The velocity-Verlet integrator46 was utilized to integrate the equations of motion. The Nosé-Hoover thermostat51-52 with a relaxation time of 0.1 ps was used to control the temperature. The Parrinello-Rahman barostat53 with a coupling time of 1.0 ps was chosen to keep the constant pressure. In fact, a new equilibrium had been achieved at about 0.1 ns during the annealing process. Finally, the relaxed size of the simulation box was 52.2 Å × 52.2 Å × 52.2 Å. After the annealing step, we performed another equilibrium simulation for 10 ns under NVT ensemble at the constant temperature of 298 K, and one trajectory of all-atom model was stored in every 0.1 ns. The 100 trajectories from the allatom model under the NVT ensemble were used for CG statistics of the structural information, including the bond length, bond angle, bond dihedral, and nonbonded RDF. Besides, the relaxed all-atom model of PI was stretched at 1010 s–1 strain rate in X direction, and the engineering strain was up to 1.2; the stress-strain relation was recorded.

Figure 2. Mapping scheme of poly(cis-1,4-isoprene) from the full atomistic model to the coarse-

grained model.54 Gray sphere, carbon atom; white sphere, hydrogen atom; green sphere, repeating unit; blue sphere, superatom; , center of the superatom. The (N – 1)-mer scheme,54 proposed by Reith at al.17, 55 for the mapping of poly(trans-1,4isoprene) and later used for the mapping of poly(cis-1,4-isoprene) by Faller,54 was used in the


14

Page 15 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


statistical analysis of the conformational distributions. Figure 2 shows the mapping details of the superatom in one chain of PI. The center of the superatom was located at the middle of the carbon-carbon single bond connecting the two successive repeat units, and each chemical repeat unit was sketched as one superatom. We selected the (N – 1)-mer scheme to statistically analyze the conformational distributions, because the choice of the superatom center being placed between atomistic monomers of PI can generate one distinguishable peak in bond length distribution. This was beneficial for the intramolecular distributions to increase the independence between the conformational distributions.54 The starting coordinates of the superatoms of the CG model were determined by coarsegraining the all-atom model of PI according to the (N – 1)-mer mapping scheme. The CG potentials were obtained via the standard IBI CG method. For each iteration, the standard IBI CG simulation was carried out for 10 ns, and the trajectory was recorded every 0.1 ns. Then these trajectories were used for the statistics of the probability distributions. Ubond, Uangle, Unonbonded, and Udihedral were gained after 12 times of iterations for convergence. The CG simulations, which were carried out by using LAMMPS, were divided into three categories: the first category was the standard IBI CG simulation for structural CG potentials, the second category was the standard IBI CG simulation in combination with the DPD thermostat (standard IBI+DPD), and the last category was the CG simulation of modifying the IBI nonbonded potential to compensate for the reduced friction by using the additional function. As for the second and the third categories, the bonded potentials derived from the standard IBI method were kept intact, while the nonbonded potentials were changed. In the first category, the size of the simulation box was 52.2 Å × 52.2 Å × 52.2 Å, which was identical with the relaxed all-atom model. Then, the standard IBI procedure was performed to obtain the structural CG potentials under the NVT


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 39

ensemble. The cutoff radius rcutoff, the time step, the temperature and other conditions were the same as the target all-atom model. The Nosé-Hoover thermostat51-52 with a relaxation time of 0.1 ps was used to control the temperature of the CG model; and the Parrinello-Rahman barostat53 with a coupling time of 1.0 ps was chosen to keep the constant pressure. In the second category, the nonbonded potential was overlaid with the coupled dissipative-random forces of the DPD thermostat, and the dissipative-random forces were controlled by the strength parameter γ. In our simulation, γ was specified by trying to match the tensile stress amplitude of the CG model with that of the all-atom model. As for the last category, the additional empirical function had been introduced into the nonbonded potential. Some parameters of G wave ( r , ∆t ) were determined by keeping the stress-strain results from the CG model in consistence with that from the all-atom model, while the other parameters were obtained from the thermal fluctuation information of the nonbonded RDF from the all-atom model. 2.4. Physical Property Calculation. To examine the effect of the additional empirical

function in eq 24 on the physical properties, such as the pressure, self-diffusion coefficient, glass-transition temperature, Hamiltonian and enthalpy, were computed after the CG force fields had been obtained. In our simulation, a self-diffusion coefficient, Dc , was calculated under NVT ensemble at a temperature of 298 K according to the Einstein equation:56-57 Dc = lim t →∞

g MSD (t ) , 6t

g MSD (t ) = ri ( t ) − ri ( 0 )

(26) 2

,

(27)

where the time-dependent g MSD (t ) is the mean-squared displacement of the individual superatoms, ri ( t ) is the position of the i-th superatom at time t, and operator denotes the average of the sum.


16

Page 17 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The glass-transition temperature, Tg , was obtained by calculating the temperature-dependent mean-squared displacement of the polymer segments according to the approach proposed by Tsige and Taylor.58 To search for Tg , the temperature interval of 20 K was selected. Before measuring the mean-squared displacement, the model was equilibrated at each temperature for 5 ns under NVT ensemble. The last trajectory in the current stage was considered as the starting trajectory for the calculation of mean-squared displacement. In the next stage, the relaxed model was further simulated for another 1 ns under NVT ensemble, and the last trajectory was taken as the ending trajectory for the calculation. The mean-squared displacement was calculated using eq 27. Thus, a series of the mean-squared displacements at different temperatures were reached. When these data were drawn in a plot, the variation of the mean-squared displacement with the temperature can be approximately described by two different straight lines, whose intersection demonstrated the corresponding Tg . The Hamiltonian and enthalpy, as the important thermodynamic properties, were calculated. The Hamiltonian H m was the sum of the total potential energy Pe and the total kinetic energy K e , i.e., H m =Pe + K e .

(28)

The enthalpy H Ω was calculated by H Ω = H m + PV .

(29)

3. RESULTS AND DISCUSSION 3.1. Standard IBI CG Potential. The CG potentials of PI bulk were obtained in light of the

process of the standard IBI method. Figure 3 shows the CG potentials including the bonded potentials (stretching potential, the bending potential, and the torsion potential) and the


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 39

nonbonded potential. From the comparison of the structure distributions, it can be seen in Figure 3 that the distributions of bond length, bond angle, RDF, and bond dihedral obtained from the all-atom calculation of using CVFF50 were all in good consistence with those of the CG calculation of using the standard IBI force fields, respectively. The results indicated that we had obtained the validated IBI CG potentials due to their capability to reproduce the target structural distributions. Additionally, the obtained results were in agreement with the IBI CG potentials based on the Monte Carlo simulations.59 The bonded potentials were defined in the analytical forms as listed in Table 1, for it was convenient to adjust their parameters for sensitive calculations. Furthermore, the nonbonded potential after the pressure correction was divided into several portions according to the distance r between the superatoms; and the fitting form and parameters are shown in Tables 1 and 2, respectively. In the following CG simulation, the nonbonded potential was modified to obtain the mechanical representability of the CG PI model, while the bonded potentials were used for keeping the conformational structure of the CG model.


18

Page 19 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

Figure 3. Probability distributions and potentials for PI. (a) Bond length distributions between

two consecutive CG superatoms and corresponding potential. Green thick line, distribution from the all-atom model; red thin line, distribution from the coarse-grained model; blue dot-dash line, final corresponding numerical potential derived from the standard IBI method. (b) As in (a) for bond angle θ. (c) As in (a) for nonbonded RDF and corresponding potential. rcore is the radius of the excluded volume of the nonbonded superatom; rfirst is the corresponding distance r of the first peak value of RDF. (d) As in (a) for bond dihedral angle φ .


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 39

Table 1. Functional Forms of Coarse-Grained Force Fields and Potential Parameters for PI

interaction

potential form

parameters

bond length

U bond ( l ) = k ( l − l0 ) 2

k = 4.49 kcal/(mol Å2) l0 = 3.863 Å

bond angle

2  1   θ − θ1   a1 exp −  U angle (θ ) = − k B T ln    b1   sin(θ ) 

 θ − θ 3     θ − θ2  exp + a 2 exp −  + a −     3  b2   b3    2

bond dihedral

2

3

Udihedral (φ ) = ∑ Ai cosi −1 (φ ) i =1

a1 = 0.3698 a 2 = 0.9680 b1 = 29.44 o b2 = 45.12 o

θ1 = 40.05 o θ 2 = 89.99 o a3 = 0.3713 b3 = 29.46o

θ3 = 139.94o A1 = 3.463 kcal/mol

A2 = −0.1956 kcal/mol A3 = 0.07299 kcal/mol

nonbonded

2  n  r − ri   U nonbonded (r ) = −k BT ln ∑ ai exp −     i =1  bi  

see Table 2

3.2. Determination of the Additional Empirical Parameters. The parameters of additional

empirical function G wave ( r , ∆t ) were obtained by abstracting the thermal fluctuations of the nonbonded RDF of the all-atom model. Subsequently, two crucial parameters in the additional empirical function were determined by adjusting the amplitude of G wave ( r , ∆t ) to make sure that, under the same stretching condition, the stress-strain curve of the CG model was in consistence with that of the all-atom model. In detail, the initial determination of f outline ( r ) was based on fitting the stochastic amplitudes of the thermal fluctuations of the time-dependent nonbonded RDFs. We at first statistically analyzed one hundred trajectories of the all-atom simulation under NVT ensemble to obtain individual time-dependent RDF of the nonbonded superatom. Subsequently, the mean of all the time-dependent RDFs was just the time-averaged RDF, shown


20

Page 21 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Figure 3c. From comparison of the one hundred time-dependent RDFs with the time-averaged 1 2 100 RDF, the set of differences Ddiffer (r ) = ( Ddiffer , Ddiffer ,K , Ddiffer ) , i.e., the thermal fluctuations of T

the nonbonded RDF, were collected by subtracting the time-averaged RDF from any of the timedependent RDFs. Table 2. Nonbonded Potential Parameters of the Standard IBI Method

distance scope r/Å

parameters n

r1 /Å

a1

b1

r2 /Å

a2

b2

r3 /Å

a3

b3

/

/

/

/

(0, 4.6]

1 5.297

0.930

0.990

/

/

(4.6, 6.3]

3 5.051

0.646

0.869

5.085

0.0687 0.319

6.893

0.801

1.443

(6.3, 7.8]

3 5.840

0.477

0.633

6.933

0.529

0.864

9.706

1.336

2.362

(7.8, 9.3]

3 7.484

0.0777 0.241

8.355

0.922

2.351

14.506 29.773 2.240

(9.3, 10.8]

3 9.847

0.0069 0.0707 7.429

0.866

1.708

10.827 0.976

2.433

(10.8, 15]

3 10.164 0.066

14.833 0.0301 1.600

21.903 0.998

41.733

2.590

Figure 4 shows that the differences Ddiffer (r ) can be stacked in an area range with a clear i outline, even though the difference at a specific moment Ddiffer (r ) with i = 1, 2, 3,...,100 has a

strong randomness that is a typical sign of the thermal fluctuations. The fluctuations in the differences

Ddiffer ( r )

Ddiffer ( r ) = Ddiffer

1

of

the

nonbonded

super-atom

were

averaged

over,

as

in

100 . And then, the data of Ddiffer ( r ) within (rfirst, rcutoff) were fitted via an

exponential function Ffit ( r ) for the original trend of the outline function f outline ( r ) .


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 39

Figure 4. Thermal fluctuations of the RDF of the nonbonded superatom from the all-atom model.

The cyan dots are the stacks of all the elements of Ddiffer (r ) , and the blue curve is the fluctuation i information of Ddiffer (r ) at i = 1.

In eq 25, η ( r ) is a probability ratio defined as,

η (r ) = g 0 (r ) / g 0 (rfirst )，

(30)

where rfirst = 5.24 Å in our simulation. f outline ( r ) in the sinusoidal function of eq 25 is an adjustable piecewise function. For rfirst ≤ r ≤ rcutoff , f outline ( r ) was given by f outline ( r )  π ( r − radd )  = C1 Ffit (r )+C 2 Ffit (rfirst ) sin  , 2 r − r ( )  cutoff add 

(31)

in which C1 and C2 are two parameters for adjusting the amplitudes of the outline function

f outline ( r ) ; radd on behalf of the beginning point of the sinusoidal function in eq 31 was set to ( rfirst + rcutoff )/2, and the value of C2 with r in the range of ( rfirst , radd ) was always fixed to be zero in any CG simulation. In the range of ( rcore , rfirst ), f outline ( r ) was given by f outline ( r ) = C3 ( r − rcore ) ( rfirst − rcore ) 6

−6

to keep the trend in consistence with the steeper outline in Figure 4;


22

Page 23 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rcore = 2.96 Å in our simulation. The constant C3 was determined by using the continuous

condition at r = rfirst . In the range of (0, rcore ), foutline ( r ) = 0 , because there was no fluctuations i (r ) in this range shown in Figure 4. Besides, in the range of ( rcutoff , + ∞ ), the for Ddiffer

contribution of the fluctuations of nonbonded RDF to the nonbonded potential was assumed to be zero in the present work. In additional empirical function G wave ( r , ∆t ) , rcycle , the mean distance length of the fluctuation waves near rfirst , was obtained from the statistical analysis of all the elements of Ddiffer (r ) ; here

rcycle = 0.165 Å in our present work. ϕ 0 was on behalf of the original phase angle of G wave ( r, ∆t ) at

rfirst . In principle, ϕ 0 should be given by a time-dependent value. For the simplicity, ϕ 0 =

3π /2 in the present simulation.

In the range of (0, rcore), the function g(r) from Gaussian fitting of the RDF of the nonbonded superatom was greater than zero, and the potential values from the Boltzmann inversion were much smaller than positive infinity. Thus, the volume-excluded core of the nonbonded superatom under the force fields was softer than the spherical superatom core from the all-atom model, which may result in excessive collapse of the volume of the CG model under fast strainrate stretching. In order to avoid this unrealistic phenomenon, we introduced a correction factor δ

( rcore / r )

into the nonbonded potential function in the range of (0, rcore) for hardening the

volume-excluded core, as follows, δ

r  DPD,core U nonbonded (r ) =  core  U nonbonded (r )  r   G wave ( r , ∆t )  − kBT ln 1 + , g IBI ( r )  

0 < r < rcore

(32)


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 39

where δ ≥ 1 , and the greater the value of δ, the harder the volume-excluded core was; here

δ = 7 .60 It was found that the modification of the excluded volume core of the nonbonded superatom had little influence on the pressure of the CG model. Thus there was no need to readjust the pressure. According to the above process, some parameters of the modified nonbonded potential can be determined in the case of the thermal fluctuation information into additional empirical function. The other two important parameters, C1 and C2, were adjustable. In order to figure out the effect of the two parameters on the mechanical property, a series of tensile tests under stretching in X direction with a strain rate of 1010 s–1 were carried out. As shown in Table 3, the maximum stress

σ max increased with an increase of any one of the two parameters. Table 3. Maximum Stresses σ max at Different Parameters

σ max (MPa) C1 C2

8.7 0.5 0

16.9 45.3 153.6 14.5 28.4 110.0 468.0 1.0 2.0 4.0 0 0 0 0 0 0 0 0.5 1.0 2.0 4.0

In view of the effect of the two parameters on the mechanical property, an approach of optimizing C1 and C2 was proposed by matching the stress-strain relationship of the all-atom simulation under uniaxial stretching. At first, let C2 be zero in this example, select the maximum stress of the all-atom model, Kpeak, as the target key point; and adjust C1 to tune the stress at strain corresponding to Kpeak of the CG model to be close to that of the all-atom model. Subsequently, C2 was gradually increased, and C1 was simultaneously decreased a little to recover the target peak of the stress-strain curve for the combination C1 and C2 will lead to an increase of the maximum stress if C1 kept intact. The obtained parameters can be evaluated by the following evaluation function,


24

Page 25 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∫ [σ (ε ) − σ (ε )] dε , = ∫ [σ (ε )] d ε 2

f eval

target

(33)

2

target

in which σ (ε ) and σ target (ε ) are the stress functions from the CG model and from the target allatom model, respectively, and ε is strain. Obviously, the smaller value of the evaluation function, the better result will be obtained for the two parameters. For example, f eval = 0.0926 in the case of C1 = 4.30 and C2 = 0; When C1 = 3.92 and C2 = 0.43, f eval = 0.0545 . So, the latter was a better result for the two parameters. 3.3. Comparison of CG Models and All-Atom Model. The mechanical property of the CG

model is an essential aspect. We compared the simulation results of the all-atom and the CG models of PI under stretching at 1010 s–1 strain rate in X direction. The tensile results show in Figure 5 that the stress-strain relation of the standard IBI model is the lowest in amplitude among the three CG models. This revealed that the standard IBI method sometimes cannot reproduce the correct mechanical property. The similar mechanical results were also found for the standard IBI CG simulation of polystyrene.26-27 Müller-Plathe et al.,20, 26 attributed the notable difference of the mechanical property mainly to the much softer (less attractive) and smoother CG potential in comparison with the all-atom potential. In a different view, it was considered as the results from the reduction of the intrinsic friction force during the coarse-graining process.9, 24


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 39

Figure 5. Stress vs. strain curve under stretching in X direction with a strain rate of 1010 s–1. The

green square

curve shows the result from the all-atom model; the cyan diamond

curve shows

the result from the CG model using the standard IBI force fields; the magenta circle

curve

represents the result from the CG model using the standard IBI potentials in combination with the DPD thermostat; the blue triangle

curve is the result from the CG model using the

modified nonbonded force field combining with the standard IBI method. The additional function was directly related to the nonbonded potential of the modified IBI model, while the standard IBI+DPD model tuned the nonbonded potential by the DPD thermostat in an indirect process, so there existed a difference in the tuning parameters between the two models. For example, C1 and C2 for the former CG model were the amplitudes of the thermal fluctuations of the nonbonded RDF while γ for the latter CG model was the strength of the random force. All the parameters can be tuned to match the mechanical properties of the allatom model. It was found in Figure 5 that the modified IBI CG model with introduction of the additional empirical function (C1 = 3.92 and C2 = 0.43) and the standard IBI+DPD CG model (γ = 12.50 kcal·ps·mol–1·Å–2) can reproduced the strain vs. stress curve of the all-atom model in a wide range of strain. The values of the evaluation function calculated from the modified IBI and the standard IBI+DPD CG models were 0.0545 and 0.1318, respectively. Even though there was a little discrepancy in the tensile stress between the proposed CG and the standard IBI+DPD models, the modified IBI CG model can replicate the effect of the DPD thermostat’s dissipativerandom forces. This implied that friction-force reduction of the standard IBI model was compensated for by introducing the additional empirical function into the standard IBI nonbonded potential. Furthermore, it can be demonstrated by the stress error ratio of the difference between the CG stress and the all-atom stress to the CG stress. The stress error ratio vs.


26

Page 27 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


strain curve in Figure 6 indicated that, to a certain extent, the developed CG force fields extended the standard IBI CG model to the mechanical applications.

Figure 6. Stress error ratio vs. strain curve under stretching in X direction with a strain rate of

1010 s–1. The cyan diamond

curve shows the result from the CG model using the standard IBI

force fields; the magenta circle

curve represents the result from the CG model using the

standard IBI potentials in combination with the DPD thermostat; the blue triangle

curve is the

result from the proposed CG model. In view of the good independence of the nonbonded potential on the bonded potential, more attention was paid to the nonbonded RDF after modification of the CG nonbonded potential.11, 61 The structural comparison of the CG and the all-atom models was used to assess whether the modification of the nonbonded potential had an effect on the reproduction of the time-averaged nonbonded RDF. The new statistical RDF of the nonbonded superatom was plotted in Figure 7, indicating that the time-averaged nonbonded RDF of the CG model was consistent with that of the all-atom model; so the modification of the nonbonded RDF by introducing the additional empirical function had little altered the local packing structure of the superatom chains of the CG model. Consequently, using the nonbonded force field modified by the additional empirical


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 39

function combining with the IBI bonded force fields, the CG model was still capable of reproducing the structural information.

Figure 7. Time-averaged RDFs of the nonbonded superatom. The green thick curve is from the

all-atom model; the red thin curve is from the CG model using modified nonbonded potential combining with the standard IBI method. The thermodynamic properties, such as the pressure, self-diffusion coefficient, and glasstransition temperature, were calculated using the all-atom model, the standard IBI CG model, the standard IBI+DPD CG model, and the modified IBI CG model under NVT ensembles, as listed in Table 4. It can be seen that the glass-transition temperature Tg from the standard IBI+DPD model and from the modified IBI CG model were 195 K and 194 K, respectively, close to the target, 200 K from the all-atom model. Clearly, the modified IBI CG model as same as the standard IBI+DPD model was better for matching the glass transition than the standard IBI model, whose Tg was 84 K. Since the modified IBI CG model had the same technique for pressure correction as the standard IBI CG model did, the pressure calculated from the former was 3.6 atm, almost equal to that from the latter, 3.2 atm; the pressure of the standard IBI+DPD model was 12.5 in our simulation, indicating that the DPD thermostat had a little effect on the pressure. As for the comparison of such thermodynamic properties as the self-diffusion


28

Page 29 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


coefficient Dc , Hamiltonian H m and enthalpy H Ω , there existed considerable discrepancies between the standard IBI CG model and the all-atom model. The introduction of the DPD thermostat into the standard IBI model had little effect on the H m and H Ω . After incorporation of the additional empirical function into the standard IBI nonbonded potential, all the discrepancies were decreased more or less. Even that, those thermodynamic properties were not well reproduced by the modified IBI CG model; the reason might be that the parameters of the additional empirical function were merely obtained from matching the mechanical properties rather than the thermodynamic properties. More efforts will be done for the improvement of the modified IBI CG model to recover more physical properties of the all-atom model. Actually, it is difficult to keep all the physical properties in simultaneous consistence between the all-atom model and the CG model.9, 28 Table 4. Thermodynamic Properties Calculated from Four Models

Model all-atom standard IBI CG standard IBI+DPD CG modified IBI CG

P (atm) 1.7 3.2 12.5 3.6

Tg (K) 200 84 195 194

Dc (Å2/ns) 0.51 397.12 4.09 9.33

Hm (kcal/mol) 24053 4599 4607 9730

HΩ (kcal/mol) 24037 4606 4809 9750

4. CONCLUSIONS

The standard IBI nonbonded potential was modified by adding an empirical function. According to an analysis of the dissipative particle dynamics, the additional function did compensate for the friction reduction of the standard IBI CG model. To demonstrate the nonbonded potential-modified method, an additional empirical function of the thermal fluctuation information from the nonbonded radial distribution function of poly(cis-1,4-isoprene) was selected as an example. The simulation results indicated that the CG force fields of the


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 39

method could ensure that the stress-strain relation and the glass-transition temperature were consistent with those from the target all-atom model. In other words, we modified the CG nonbonded potential derived from the standard IBI method to take into account the reduced dissipative forces between the nearby nonbonded superatoms of the CG model. Additionally, there existed the discrepancies of such thermodynamic properties as self-diffusion coefficient, Hamiltonian and enthalpy between the modified IBI CG model and the all-atom model. More efforts will be done for the improvement of the modified IBI CG model to recover more physical properties of the all-atom model. AUTHOR INFORMATION Corresponding Author

* Z. Xie. E-mail: [email protected]. Phone: +86 0451 86402378. Author Contributions

†Z. Xie and D. Chai contributed equally to this work.

ACKNOWLEDGMENT The authors gratefully acknowledge financial support by the National Science Foundation of P. R. China under Projects No.10972068. The authors expressed special thanks to the reviewers’ valuable suggestions and comments. REFERENCES (1) Padding, J. T.; Briels, W. J. Systematic Coarse-Graining of the Dynamics of Entangled Polymer Melts: the Road from Chemistry to Rheology. J. Phys.: Condens. Matter 2011, 23, 233101.


30

Page 31 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2) Harmandaris, V. A.; Kremer, K. Dynamics of Polystyrene Melts through Hierarchical Multiscale Simulations. Macromolecules 2009, 42, 791-802. (3) Snow, C. D.; Nguyen, H.; Pande, V. S.; Gruebele, M. Absolute Comparison of Simulated and Experimental Protein-Folding Dynamics. Nature 2002, 420, 102-106. (4) Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; et al. Atomic-Level Characterization of the Structural Dynamics of Proteins. Science 2010, 330, 341-346. (5) Hyeon, C.; Thirumalai, D. Capturing the Essence of Folding and Functions of Biomolecules Using Coarse-Grained Models. Nat. Commun. 2011, 2, 487. (6) Irannejad, R.; Tomshine, J. C.; Tomshine, J. R.; Chevalier, M.; Mahoney, J. P.; Steyaert, J.; Rasmussen, S. G. F.; Sunahara, R. K.; El-Samad, H.; Huang, B.; et al. Conformational Biosensors Reveal GPCR Signalling from Endosomes. Nature 2013, 495, 534-538. (7) Simons, K.; Toomre, D. Lipid Rafts and Signal Transduction. Nat. Rev. Mol. Cell Biol. 2000, 1, 31-39.

(8) Peter, C.; Kremer, K. Multiscale Simulation of Soft Matter Systems - from the Atomistic to the Coarse-Grained Level and Back. Soft Matter 2009, 5, 4357-4366. (9) Hsu, D. D.; Xia, W.; Arturo, S. G.; Keten, S. Thermomechanically Consistent and Temperature Transferable Coarse-Graining of Atactic Polystyrene. Macromolecules 2015, 48, 3057-3068.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 39

(10) Ingolfsson, H. I.; Lopez, C. A.; Uusitalo, J. J.; de Jong, D. H.; Gopal, S. M.; Periole, X.; Marrink, S. J. The Power of Coarse Graining in Biomolecular Simulations. Wiley Interdiscip.

Rev.: Comput. Mol. Sci. 2014, 4, 225-248. (11) Karimi-Varzaneh, H. A.; Müller-Plathe, F. In Multiscale Molecular Methods in Applied

Chemistry; Kirchner, B.; Vrabec, J., Eds. Springer: Berlin, 2012; Vol. 307, pp 295-321. (12) Klein, M. L.; Shinoda, W. Large-Scale Molecular Dynamics Simulations of SelfAssembling Systems. Science 2008, 321, 798-800. (13) Frembgen-Kesner, T.; Andrews, C. T.; Li, S.; Ngo, N. A.; Shubert, S. A.; Jain, A.; Olayiwola, O. J.; Weishaar, M. R.; Elcock, A. H. Parametrization of Backbone Flexibility in a Coarse-Grained Force Field for Proteins (COFFDROP) Derived from All-Atom Explicit-Solvent Molecular Dynamics Simulations of All Possible Two-Residue Peptides. J. Chem. Theory

Comput. 2015, 11, 2341-2354. (14) Fritz, D.; Koschke, K.; Harmandaris, V. A.; van der Vegt, N. F. A.; Kremer, K. Multiscale Modeling of Soft Matter: Scaling of Dynamics. Phys. Chem. Chem. Phys. 2011, 13, 1041210420. (15) Brini, E.; Algaer, E. A.; Ganguly, P.; Li, C. L.; Rodriguez-Ropero, F.; van der Vegt, N. F. A. Systematic Coarse-Graining Methods for Soft Matter Simulations - A Review. Soft Matter 2013, 9, 2108-2119.

(16) Müller-Plathe, F. Coarse-Graining in Polymer Simulation: From the Atomistic to the Mesoscopic Scale and Back. ChemPhysChem 2002, 3, 754-769.


32

Page 33 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24, 1624-1636. (18) Milano, G.; Müller-Plathe, F. Mapping Atomistic Simulations to Mesoscopic Models: A Systematic Coarse-Graining Procedure for Vinyl Polymer Chains. J. Phys. Chem. B 2005, 109, 18609-18619. (19) Qian, H.-J.; Carbone, P.; Chen, X.; Karimi-Varzaneh, H. A.; Liew, C. C.; Müller-Plathe, F. Temperature-Transferable Coarse-Grained Potentials for Ethylbenzene, Polystyrene, and Their Mixtures. Macromolecules 2008, 41, 9919-9929. (20) Carbone, P.; Varzaneh, H. A. K.; Chen, X.; Müller-Plathe, F. Transferability of CoarseGrained Force Fields: The Polymer Case. J. Chem. Phys. 2008, 128, 064904. (21) Henderson, R. L. A Uniqueness Theorem for Fluid Pair Correlation Functions. Phys. Lett.

A 1974, 49, 197-198. (22) Ganguly, P.; van der Vegt, N. F. A. Representability and Transferability of Kirkwood– Buff Iterative Boltzmann Inversion Models for Multicomponent Aqueous Systems. J. Chem.

Theory Comput. 2013, 9, 5247-5256. (23) Moore, T. C.; Iacovella, C. R.; McCabe, C. Derivation of Coarse-Grained Potentials via Multistate Iterative Boltzmann Inversion. J. Chem. Phys. 2014, 140, 224104. (24) Qian, H.-J.; Liew, C. C.; Müller-Plathe, F. Effective Control of the Transport Coefficients of A Coarse-Grained Liquid and Polymer Models Using the Dissipative Particle Dynamics and Lowe-Andersen Equations of Motion. Phys. Chem. Chem. Phys. 2009, 11, 1962-1969.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 39

(25) Li, Y.; Abberton, B. C.; Kroeger, M.; Liu, W. K. Challenges in Multiscale Modeling of Polymer Dynamics. Polymers 2013, 5, 751-832. (26) Rahimi, M.; Iriarte-Carretero, I.; Ghanbari, A.; Böhm, M. C.; Müller-Plathe, F. Mechanical Behavior and Interphase Structure in A Silica–Polystyrene Nanocomposite under Uniaxial Deformation. Nanotechnology 2012, 23, 305702. (27) Rosch, T. W.; Brennan, J. K.; Izvekov, S.; Andzelm, J. W. Exploring the Ability of a Multiscale Coarse-Grained Potential to Describe the Stress-Strain Response of Glassy Polystyrene. Phys. Rev. E 2013, 87, 042606. (28) Gao, P.; Guo, H. Developing Coarse-Grained Potentials for the Prediction of MultiProperties of trans-1,4-Polybutadiene Melt. Polymer 2015, 69, 25-38. (29) Harmandaris, V. A.; Reith, D.; van der Vegt, N. F. A.; Kremer, K. Comparison Between Coarse-Graining Models for Polymer Systems: Two Mapping Schemes for Polystyrene.

Macromol. Chem. Phys. 2007, 208, 2109-2120. (30) Harmandaris, V. A.; Kremer, K. Predicting Polymer Dynamics at Multiple Length and Time Scales. Soft Matter 2009, 5, 3920-3926. (31) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K. Hierarchical Modeling of Polystyrene: From Atomistic to Coarse-Grained Simulations. Macromolecules 2006, 39, 6708-6719.

(32) León, S.; van der Vegt, N.; Delle Site, L.; Kremer, K. Bisphenol A Polycarbonate: Entanglement Analysis from Coarse-Grained MD Simulations. Macromolecules 2005, 38, 80788092.


34

Page 35 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(33) Depa, P. K.; Maranas, J. K. Speed up of Dynamic Observables in Coarse-Grained Molecular-Dynamics Simulations of Unentangled Polymers. J. Chem. Phys. 2005, 123, 094901. (34) Depa, P. K.; Maranas, J. K. Dynamic Evolution in Coarse-Grained Molecular Dynamics Simulations of Polyethylene Melts. J. Chem. Phys. 2007, 126, 054903. (35) Chen, C.; Depa, P.; Sakai, V. G.; Maranas, J. K.; Lynn, J. W.; Peral, I.; Copley, J. R. D. A Comparison of United Atom, Explicit Atom, and Coarse-Grained Simulation Models for Poly(ethylene oxide). J. Chem. Phys. 2006, 124, 234901. (36) Li, Y.; Tang, S.; Abberton, B. C.; Kröger, M.; Burkhart, C.; Jiang, B.; Papakonstantopoulos, G. J.; Poldneff, M.; Liu, W. K. A Predictive Multiscale Computational Framework for Viscoelastic Properties of Linear Polymers. Polymer 2012, 53, 5935-5952. (37) Strauch, T.; Yelash, L.; Paul, W. A Coarse-Graining Procedure for Polymer Melts Applied to 1,4-Polybutadiene. Phys. Chem. Chem. Phys. 2009, 11, 1942-1948. (38) Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A.; Kremer, K.; Mann, B. A.; Voelkel, R.; Weiss, H.; Liew, C. Ethylbenzene Diffusion in Polystyrene: United Atom Atomistic/Coarse Grained Simulations and Experiments. Macromolecules 2007, 40, 7026-7035. (39) Hsu, D. D.; Xia, W.; Arturo, S. G.; Keten, S. Systematic Method for Thermomechanically Consistent Coarse-Graining: A Universal Model for Methacrylate-Based Polymers. J. Chem.

Theory Comput. 2014, 10, 2514-2527. (40) Xia, W.; Mishra, S.; Keten, S. Substrate vs. Free Surface: Competing Effects on the Glass Transition of Polymer Thin Films. Polymer 2013, 54, 5942-5951.


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 39

(41) Ganguly, P.; Mukherji, D.; Junghans, C.; van der Vegt, N. F. A. Kirkwood–Buff CoarseGrained Force Fields for Aqueous Solutions. J. Chem. Theory Comput. 2012, 8, 1802-1807. (42) Moore, T. C.; Iacovella, C. R.; McCabe, C. In Foundations of Molecular Modeling and

Simulation: Select Papers from FOMMS 2015; Snurr, Q. R.; Adjiman, S. C.; Kofke, A. D., Eds. Springer Singapore: Singapore, 2016; pp 37-52. (43) Español, P.; Warren, P. Statistical Mechanics of Dissipative Particle Dynamics. Europhys.

Lett. 1995, 30, 191-196. (44) Groot, R. D.; Warren, P. B. Dissipative Particle Dynamics: Bridging the Gap between Atomistic and Mesoscopic Simulation. J. Chem. Phys. 1997, 107, 4423-4435. (45) Wang, H.; Junghans, C.; Kremer, K. Comparative Atomistic and Coarse-Grained Study of Water: What do We Lose by Coarse-Graining? Eur. Phys. J. E 2009, 28, 221-229. (46) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76, 637-649. (47) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155-160. (48) Nikunen, P.; Karttunen, M.; Vattulainen, I. How Would You Integrate the Equations of Motion in Dissipative Particle Dynamics Simulations? Comput. Phys. Commun. 2003, 153, 407423. (49) Kubo, R. The Fluctuation-Dissipation Theorem. Rep. Prog. Phys. 1966, 29, 255-284.


36

Page 37 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Structure and Energetics of Ligand Binding to Proteins: Escherichia Coli Dihydrofolate Reductase-Trimethoprim, A Drug-Receptor System. Proteins: Struct., Funct., Bioinf. 1988, 4, 31-47. (51) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol.

Phys. 1984, 52, 255-268. (52) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695-1697.

(53) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190. (54) Faller, R. Automatic Coarse Graining of Polymers. Polymer 2004, 45, 3869-3876. (55) Faller, R.; Reith, D. Properties of Poly(isoprene): Model Building in the Melt and in Solution. Macromolecules 2003, 36, 5406-5414. (56) Charati, S. G.; Stern, S. A. Diffusion of Gases in Silicone Polymers: Molecular Dynamics Simulations. Macromolecules 1998, 31, 5529-5535. (57) Mattozzi, A.; Hedenqvist, M. S.; Gedde, U. W. Diffusivity of n-Hexane in Poly(ethylene-

stat-octene)s Assessed by Molecular Dynamics Simulation. Polymer 2007, 48, 5174-5180. (58) Tsige, M.; Taylor, P. L. Simulation Study of the Glass Transition Temperature in Poly(methyl methacrylate). Phys. Rev. E 2002, 65, 021805.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 39

(59) Pandey, Y. N.; Brayton, A.; Burkhart, C.; Papakonstantopoulos, G. J.; Doxastakis, M. Multiscale Modeling of Polyisoprene on Graphite. J. Chem. Phys. 2014, 140, 054908. (60) Karimi-Varzaneh, H. A.; van der Vegt, N. F. A.; Müller-Plathe, F.; Carbone, P. How Good Are Coarse-Grained Polymer Models? A Comparison for Atactic Polystyrene.

ChemPhysChem 2012, 13, 3428-3439. (61) Peter, C.; Kremer, K. Multiscale Simulation of Soft Matter Systems. Faraday Discuss. 2010, 144, 9-24.


38

Page 39 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC


39

Directly Modifying the Nonbonded Potential Based on the Standard

Recommend Documents