Effect of Molecular Weight on Reactive Molecular Dynamics

Subscriber access provided by Stockholm University Library

Article

Effect of Molecular Weight on Reactive Molecular Dynamics Simulations: Guiding the Theoretical Study for Macromolecules at an Atomic Level Wu Chen, Haitao Duan, Meng Hua, Yaling Xiang, Song Chen, Sheng-peng Zhan, Kali Gu, Hongfei Shang, and Jian Li J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b06016 • Publication Date (Web): 20 Nov 2015 Downloaded from http://pubs.acs.org on November 24, 2015

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 C 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 28

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

Figure 1. Initial molecular model of UHMWPE in O2 condition: (a) UHMWPE-20 + 16 O2 molecules; (b) UHMWPE-40 + 32 O2 molecules; (c) UHMWPE-80 + 64 O2 molecules.

Figure 2. Diffusion coefficients of O2 molecules in UHMWPE with different molecular weight of atomic model.

ACS Paragon Plus Environment


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 28

Figure 3. Molecular architecture and the products distribution of oxidative reaction of UHMWPE-20 at different RMD time: (a) 275.2 ps (the time of complete consumption of O2 molecules); (b) 1000ps (the last frame).



Page 3 of 28

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


Figure 5. Molecular architecture and the products distribution of oxidative reaction of UHMWPE-40 at different RMD time: (a) 256.3 ps (the time of complete consumption of O2 molecules); (b) 1000ps (the last frame) and the 3-D view of UHMWPE-80 after 1573 K−1000 ps RMD simulation.

Figure 6. The time of complete consumption of O2 molecules and the difference of thermodynamics energy.



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

Figure 7. Variation of the number of hydrogen bonds during the RMD processes.

Figure 8. Variation of the number of –OH groups during the RMD processes.


Page 4 of 28

Page 5 of 28

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


Figure 9. Variation of the number of molecular fragments during the RMD processes: (a) RMD simulation result of UHMWPE-20; (b) RMD simulation result of UHMWPE-40; (c) RMD simulation result of UHMWPE-80.



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 28

Figure 10.Variation of the number of groups during the RMD processes: (a) RMD simulation result of UHMWPE-20; (b) RMD simulation result of UHMWPE-40; (c) RMD simulation result of UHMWPE-80.

0.2 ps

80 ps

92 ps

94ps

120 ps

240 ps

600 ps

1000 ps

Figure 11. 3-D views of oxidative reaction mechanism of UHMWPE-20 from RMD simulation.


Page 7 of 28

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


Figure 12. Scheme of reactions



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 28

Effect of Molecular Weight on Reactive Molecular Dynamics Simulations: Guiding the Theoretical Study for Macromolecules at an Atomic Level Wu Chen,†,‡ Hai-tao Duan,†,‡ Meng Hua ,§Ya-ling Xiang,†,‡ Song Chen,†,‡ Sheng-peng Zhan,†,‡ Ka-li Gu, †,‡ Hong-fei Shang∥ and Jian Li*,†,‡ †

Wuhan Research Institute of Materials Protection, 126 Bao Feng Erlu, Wuhan, Hubei 430030, China

‡

State Key Laboratory of Special Surface Protection Materials and Application Technology, 126 Bao Feng Erlu, Wuhan,

Hubei 430030, China §

MBE Department, City University of Hong Kong, 83 Tat Chee Avenue, Kowloon 999077, Hong Kong

∥

State Key Laboratory of Tribology, 9003 Building, Tsinghua University, Haidian District, Beijing 100084, China

ABSTRACT: Some constructions of molecular model for reactive molecular dynamics (RMD) simulations have neglected the effect of molecular weight on the RMD results. A computational methodology for generating structural models of amorphous ultra-high molecular weight polyethylene ( UHMWPE) with different degree of polymerization was thus proposed. The methodology was then used in RMD simulations of three types of molecular models for UHMWPE, respectively with 20, 40 and 80 vinyl monomers. The accuracy of the model has been discussed by comparing the RMD results with experimental data previous ly available in literature. Results have confirmed that the atomic model with larger molecular weight allows (i) the exact calculations of diffusion coefficient and (ii) the accurate prediction of cross-linked reaction. The RMD simulations provide the basic data for understanding and for further studying chemical reaction of the other interesting macromolecules.

1. INTRODUCTION Integrative investigation of the dynamics of macromolecules by combining computer simulation with experimental study has become a research trend in recent years. Simulation of reactive molecular dynamics (RMD) is a crucial technique to study the chemical reactions in atomic level. The technique facilitates the derivation of a series of important information like molecular configuration, atomic velocities, and cleavage and formation of chemical bonds, etc. Results of RMD provide beneficial view for deeper understanding of the process and mechanism of chemical reaction for the macromolecules, which may not easily be observable in experiment due to their instantaneous and elusive nature. However, there are some limitations of RMD simulations when the condition of simulation is compared 1


Page 9 of 28

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


with its experimental counterpart. Typical examples can be related to the constraint of (i) the simulation temperature and (ii) the molecular weight of atomic model. Aiming at overcoming the imposed constraint of computational simulation time scales, many researchers incline to choose much higher temperature in their simulation rather than using the experimental temperature when they are 1-5

investigating reaction systems.

However, the use of higher temperature may affect the reaction pathways as it may

promote the reactions with high energy barriers.1 The feasibility of such simulation with higher temperature has been proved by the simulation and experimental study of Salmon et al.6 were the members of a leading team in developing a reactive force field called ReaxFF, and applied ReaxFF to study the effect of simulation and experimental temperature difference on the thermal decomposition of brown coal. Results of their simulation indicated that, despite the higher applied temperature, RMD results can still provide a good connection to experimental data within the range of temperatures. Additional to the influence of simulation temperature, molecular weight is considered as another main factor to measure the important performance index of macromolecule. Although the effect of this factor has not yet been fully explored, many researchers

7-13

use a small segment model instead of macromolecular chain to build atomic model for

RMD simulation, which is deployed specifically for investigating the decomposition and oxidation process of macromolecules. As majority of their results are mainly empirical in nature and supplemented with many experimental data, the effect of molecular weight of atomic model on RMD is still largely uncertain. Aiming at exploring theoretically the effect of molecular weight of atomic model on the data of RMD simulation, this study used ReaxFF14 RMD to simulate the chemical reaction of the macromolecule of Ultra-high molecular weight Polyethylene (UHMWPE) in the high temperature of about 1573K and oxygen-enriched environment. All the construction of atomic model and the RMD simulations were carried out in Materials Studio software.15 To understand whether the properties of UHMWPE are sensitively affected by the change of molecular weight, the RMD simulation was performed with simulation time to be as short as 1000 ps. Series of RMD simulations were carried out using UHMWPE atomic models with different molecular weight. Simulation results serve to provide theoretical guidance and basic data in the development of RMD simulation models for the research in the field of polymers. They also pave way to generate new research ideas so as to shed light for understanding the chemical reaction process of such macromolecules. 2. DETAILS OF COMPUTATION METHOD In this study, a molecular model of single-chain terminated with –CH3 groups and having the number of 20, 40 and 2



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 28

80 vinyl monomers, respectively, for UHMWPE-20, UHMWPE-40 and UHMWPE-80 was used in the simulations. 3

The density of 3-D periodic amorphous cell of the polymer chains was suitably set at 0.96g/cm . A simulation technique for obtain a global minimum energy of a system, 16 which mainly involves the exploration of “annealing dynamics” with suitably canonical ensemble of NVT (keeping constant for the number of particles (N), volume (V) and temperature (T)), was used to equilibrate the initial structures in simulations. The technique uses the force field of Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) 17 in the “annealing dynamics” process with “annealing” temperature cycle in range of 298 K and 600 K, which is controlled to fluctuate within ±10 K by using direct velocity scaling algorithm on the basis of eq 1 below.

v i ,new 

T0 v i ,old Tins tan

(1)

where vi,new and vi,old are the velocity vectors of atom I corresponding to before and after adjusting temperature, T0 is the target temperature and Tinstan the instantaneous temperature. Although direct velocity scaling can suppress the natural fluctuations of a system, it cannot be used to generate realistic thermodynamic ensembles. However, it can be controlled to accomplish quickly the equilibrium of a system, basically as an appropriate thermostat in controlling the thermaodynamic situation. Simulation in this study used (i) atom-based summation method18 to predict the interactions of Van der Waals with a cut-off distance of 16.5 Å, and (ii) Ewald method19 to reveal the interactions of long-range electrostatics between atoms. It was then followed by packing the “annealed” structure with abundant O2 molecules. To achieve this, a generalized gradient approximation (GGA-BLYP) 20 was used to optimize the O2 molecule and the parameter of spin multiplic ity to be calculated was preliminarily set to triplet. Subsequently, the following Maxwell-Boltzmann equation was applied to relate the temperature in the system to the distribution of atomic velocities. 3



mv 2 

 m  2   2 kBT  f (v)dv   4 v 2dv  e  2 k BT 

(2)

The afore well-known formula effectively correlates the probability density f(v)dv of an atom to its mass m having velocity between v and v+dv in a system at temperature T, and also to the Boltzmann constant k B . The ball-and-stick model of the initial atomic and bonding structures in UHMWPE with 10.07 Å, 12.58 Å and 15.79 Å periodic cubes, respectively, is illustrated in Figure 1. Obviously, the initially structural models consist of (i) carbon atoms as shown in green, (ii) oxygen atoms as shown in red, and (iii) hydrogen atoms as shown in grey.

3


Page 11 of 28

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


Figure 1. Initial molecular model of UHMWPE in O 2 condition: (a) UHMWPE-20 + 16 O 2 molecules; (b) UHMWPE-40 + 32 O 2 molecules; (c) UHMWPE-80 + 64 O 2 molecules.

With such models, the dynamically oxidative reaction could effectively be simulated into two distinctive stages: (i) firstly taking MD simulation time as 100 ps for preliminarily equilibrating the kinetic and potential energy distributions; and (ii) undertaking RMD simulations with canonical (NVT) ensemble and using the potential and veloc ity distribution as obtained from the equilibration stage in (i) above. The temperature was controlled using a Nosé-Hoover thermostat

21

with a Q ratio of 0.01. The main idea behind Nosé-Hoover dynamics is that a fictitious degree of freedom is added to the structure, to represent the interaction of the structure with the heat bath. The dynamic variation of potential and kinetic energy was evaluated using ReaxFF under 1000 ps RMD duration and an integral time step of 0.2 fs. ReaxFF is a precise reactive force field which was parameterized through optimization against an extensive training set obtained from Quantum mechanical (QM) calculations. The main parameters for RMD simulations of Hydrocarbon oxidation in ReaxFF included the following data: (i) dissociation energies for various bonds containing carbon, oxygen, and hydrogen, (ii) valence angle parameters for all C/H/O molecules, (iii) rotational barriers for such C-C-C-O dihedral angles, (iv) charge distributions (ReaxFF), and (v) a number of reaction energies and mechanisms relevant to hydrocarbon oxidation. 3. RMD RESULTS AND DISCUSSION 3.1. Diffusivity of O2 in UHMWPE. The ability of small molecules penetrating into polymers can be assessed by means of comparing their corresponding diffusion coefficient. The diffusivity of O2 usually varies with the molecular architecture, system temperature, and reactive time. Presence of oxidation in polymeric components is always 4



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 28

unavoidable mainly due to the attribution of diffusing oxygen into UHMWPE in the irradiation process, however, oxidative reaction (R•+O2 ) is normally controllable by suitably regulating the oxygen diffusion rate.

22

Due to (i)

almost the same products are formed in both post-irradiation and thermal oxidation processes and (ii) the level of their relative abundance is also approximately falling in a same range,

23

the diffusion coefficient of O2 in UHMWPE is thus

a major parameter dominating the proceeding characteristics of thermal-oxidation process. The use of correct diffusion coefficient of small molecules in polymer is noteworthy to result in accurate RMD simulation. 24

According to Jost and Biswas et al, such diffusion coefficient D can be mathematically expressed as eq 3 below by considering the increase in mean-square displacement (MSD) of atoms as a function of time:

D

2 1 d N lim  ri  t   ri  0  6 N t  dt i 1

(3)

where N is the number of diffusive oxygen atoms in the system; ri(t) is the positional vector of ith particle at the instant of time t; ri (0) is the originally positional vector of ith particle at the time of starting. The right hand side expression is the mean square displacement, which is described by an averaging operator ˂•˃. The diffusion coefficient D for any individual diffusive atom when the time t is sufficiently large can simply be expressed as:

D

r t   r  0 6t

2



SMSD 6t

(4)

where SMSD is a curve fitted through all MSD data against time (Figure 2).

Figure 2. Diffusion coefficients of O 2 molecules in UHMWPE with different molecular weight of atomic model.

5


Page 13 of 28

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


According to the definition in eq 4, diffusion coefficient D can be treated as one sixth of the gradient of curve of MSD data versus time. Hence, applying linear regression approach for the data of MSD and time allows the obtainment of a linear function y = kx+b which gives a best fit to all the corresponding data points, as those green straight lines in 2

2

Figure 2. The slope k of the straight lines has unit in Å /ps or cm /s, which represents the S MSD /t in eq 4. Subsequently, the diffusion coefficient D can be simplified as:

D 

k 6

(5)

The afore-illustration suggests that the diffusion coefficient D can be obtained by directly solving for the function of SMSD . In this study, the diffusion coefficient D of O2 in UHMWPE-20, UHMWPE-40 and UHMWPE-80, respectively, which was simulated with 50 ps MD duration under a room temperature of 298 K was calculated. Calculations in this study were generally undertaken under the ensemble of constant-volume and constant-energy (NVE). This was because, it varied with which thermostat method was selected to use. Furthermore, NVE dynamics normally did not artificially interfere with the thermodynamics of the system. This study recorded instantaneously the full information like temperature, energies, velocities, etc. specifically for the calculation of mean square displacement. Figure 2 displays the curve of SMSD and the diffusion coefficient D, respectively, of O2 in UHMWPE. Larger D usually means greater capacity of O2 to penetrate into UHMWPE bulk. For acquiring near realistic simulation results, the average value of diffusion coefficient D was obtained from each of the three separate 50ps MD runs. The so simulated diffusion coefficients (D20 , D40 and D80 of O2 ) for UHMWPE-20, UHMWPE-40 and UHMWPE-80, respectively, and the experimental DI and DII as reported by Daly25 and Pauly, 26 correspondingly, were listed in Table 1. Obviously, the relatively computational orientated values D40 and D80 seemed closely comparable with the average value of 1.07×10-7 cm2/s approximately for DI and DII. However, the D20 with averaged value of 2.20×10-7 cm2 /s seemed to deviate widely from the two experimental values DI and DII, and from the other two computationally simulated values D40 and D80. The large deviation of diffusion coefficient D20 from the other Ds might be due to the inaccurate atomic model with small molecular weight, typically resulting in the incorrect simulation of oxidative reaction (R•+O2 ) or the faulty analysis of molecular architecture, etc. Hence, successful RMD simulation of polymers usually involves the use of atomic model with larger molecular weight in the aspect of simulated calculation of diffusion coefficient.

6



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 28

Table 1. Diffusion Coefficient of O2 in UHMWPE Model Model

Symbol

Value (cm2 /s)

Reference

O2 in UHMWPE-20

D20

(2.20±0.09)×10-7 , Averaged

MD Calculated

-7

O2 in UHMWPE-40

D40

(1.04±0.03)×10 , Averaged

MD Calculated

O2 in UHMWPE-80

D80

(0.98±0.02)×10-7 , Averaged

MD Calculated

-7

Report I

DI

1.14×10

Ref.25

Report II

DII

1×10-7

Ref.26

3.2. Variation of molecular architecture. The molecular architecture of UHMWPE-20 and the distribution of its oxidation products are shown in Figure 3. Such structure and distribution were observed in the RMD simulation of UHMWPE-20 in an oxygen-enriched environment having temperature of 1573 K. The simulation was at a time of 275.2 ps when O2 molecules were completely consumed. The molecular architecture (Figure 3a) of such simulation suggested that: (i) there were about 10 hydrogen atoms to be oxidized by O2 and subsequently replaced by hydroxyl groups which were then connected to the carbon atoms in UHMWPE molecule; (ii) there existed with several small molecular chain fragments with backbone of C2 (i.e. two carbon atoms) or C5 (i.e. five carbon atoms) in oxidation products; and (iii) high susceptibility to have chain scissions between two carbon sites for the introduction of hydroxyl groups. The introduction of hydroxyl group subsequently resulted in destabilizing the C−C single bond. 27 Aiming at exploring the possible involvement of any further reaction, the last frame of molecular architecture of UHMWPE-20 at the instant of 1000ps in the simulation was demonstrated in Figure 3b. The architecture showed (i) no sight of any small molecular fragment in the products, and (ii) sight of forming larger cross-linked molecule with the decrease in hydroxyl groups. Such phenomena allowed speculating that high-temperature environment would activate the oxidative and cross-linked reaction process of UHMWPE in two ways: (i) thermal and oxidative decomposition to introduce hydroxyl groups and to form small molecular fragments; and (ii) formation of reticular macromolecules from chemical cross-linking for graft-polymerization.

28

These high temperature reaction processes were also deduced from

the simulated RMD results of UHMWPE-40 (Figure 4) and UHMWPE-80 (Figure 5). Hence, common process for surface treatment of UHMWPE materials often used the process of oxygen plasma to obtain cross-linked surface. 29-31

7


Page 15 of 28

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 corresponding molecular architecture of UHMWPE-40 (Figure 4) and UHMWPE-80 (Figure 5) simulated in an environment with temperature of 1573 K and RMD at time of 1000 ps appeared to be similar to the result of UHMWPE-20. Analys is of the architecture in Figures 4a and 5a facilitated the identification of (i) several small C3~C7 molecular fragments in the oxidation products and (ii) a significant increment in hydroxyl groups. In addition, some other decomposed products like formic acid and carbon dioxide were also formed and observed in the oxidation process of UHMWPE-40 and UHMWPE-80. A few carbocycles like three-, five- and six-membered rings (Figures 4b, 5a and b) were also generated during the oxidative reaction of UHMWPE-40 and UHMWPE-80. The unstable ring structures of these carbocycles were susceptible to further oxidization and branched to C-C bond which was broken from the introduction of oxidative functional group like hydroxyl. The occurrence of such unstable ring structures was hardly observed in the RMD simulation of UHMWPE-20. It was thus anticipated that the aforementioned carbocycles, as the transitional structures, would have higher chance to take place in RMD with atomic model of larger molecular weight.

8



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 28


Comparing the original UHMWPE structure (Figure 1c) with the 3-D view of UHMWPE-80 simulated in the environment with temperature of 1573 K and RMD at the instant of 1000 ps (Figure 5b), it showed that the simulation conditions of RMD inclined to randomize severely the originally well-organized structure of molecular chain. Such severe randomization was mainly due to the nature of (i) interactions of intermolecular hydrogen bonds and (ii) “re-polymerization” process, which bonds chain fragment of the terminating carbon atom (see: Bonding sites I and II, respectively, in Figure 5b) with corresponding chain fragment belonging to adjacently periodic cells . The 3-D view of typical “re-polymerized” backbones is illustrated by the two periodic amorphous cells and displayed as ball-and-stick model as shown in Figure 5b. The result seemed agreeable with our earlier RMD simulations for polymers.

32

Study indicated that a few small C5~C7 molecular fragments were still remaining in oxidation products of UHMWPE-40 and UHMWPE-80 at the instant of 1000 ps, suggesting the simulation time longer than 1000 ps to be needed for completing the cross-link fully in the RMD simulation of UHMWPE having 40 monomers or more in an atomic model.

9


Page 17 of 28

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


Figure 5. Molecular architecture and the products distribution of oxidative reaction of UHMWPE-40 at different RMD time: (a) 256.3 ps (the time of complete consumption of O2 molecules); (b) 1000ps (the last frame) and the 3-D view of UHMWPE-80 after 1573 K−1000 ps RMD simulation.

3.3. Quantitative analysis of RMD results. The change of thermodynamic state in molecular model can be assessed by the difference of total energy of molecular systems before and after RMD simulation. Calculation of the difference of total energy was carried out with the average total energy estimated for every three separate RMD runs undertaken for each type of UHMWPE molecular model. Subsequently, the difference of total energy was calculated by the following formula:

E total  E RMD  Eoriginal

(6)

where ΔEtotal is the difference of total energy of molecular system before and after RMD simulations, ERMD is the total 10



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 28

energy of molecular system after RMD simulations, and Eoriginal is the total energy of the initially molecular system. The plot for the estimated values (Figure 6) suggested that doubling molecular weight of atomic model doubles the difference of total energy, although the RMD s imulation time for completely consuming all the O2 molecules for UHMWPE-20, UHMWPE-40 and UHMWPE-80 were almost the same, with value averagely about 275 ps. Hence, effect of molecular weight of atomic model on (i) thermodynamic performance and (ii) chemical reaction rate, respectively, was negligible in the RMD simulations for such homopolymers.

Figure 6. The time of complete consumption of O 2 molecules and the difference of thermodynamics energy.

The number of hydrogen bonds and –OH groups in the simulations varying with RMD simulation time is shown in Figure 7 and Figure 8, respectively. As there was not any hydrogen bond in the initial state of the molecular systems (Figure 1), hydrogen bond was formed and its number increased with the formation of (i) hydroxyl groups and (ii) H2 O molecules when RMD simulation time was proceeding (Figure 5). The changing number of hydroxyl groups (Figure 8) during the RMD simulations for the three types of UHMWPE molecular models tended to resume similar trends with RDM simulation time, typically as: (i) rapidly increasing the number of –OH groups with RMD simulation time shorter than 220 ps approximately; and (ii) slowly decreasing with RMD simulation time at or longer than 220 ps, which subsequently resumes a trend to stabilize at the last few picoseconds of the individual simulations. Such changing trend of the number of –OH groups (Figure 8) indicated that the formation of hydrogen bonds was mainly attributed to the interactions between hydroxyls or between hydroxyls and water molecules before 220 ps. Although there was still increasing in the number of hydrogen bonds (Figure 7) after 220 ps, the number of –OH groups was gradually decreased to a stabilization stage in the last few picoseconds. The gradual decrease in –OH groups after 220 11


Page 19 of 28

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


ps implied that the subsequent formation of hydrogen bonds in the RMD simulation period from 220 ps to 1000 ps was mainly dominated by the increase in the number of water molecules.

Figure 7. Variation of the number of hydrogen bonds during the RMD processes.

Figure 8. Variation of the number of –OH groups during the RMD processes.

Figure 9 shows the number change in the main products during the RMD simulation processes. Common trends among the quantitative characterization of oxidation products in the three types of UHMWPE models (Figure 9a, b and c) were observed in the environmental temperature of 1573 K and RMD simulation time of 1000 ps. They are typically: (i) with increasing time, the number of oxidation products of small C3~C4 and C5~C7 molecular chains increased firstly and subsequently, from about 200 ps onwards, tended to decrease to a stabilization stage; and (ii) the number of H2 O molecules tended to be doubled when the molecular weight of atomic model was doubled. Analys is suggested negligible effect of molecular weight of molecular model on the process of oxidative reaction in dehydrogenation, 12



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 28

formation of C=C bond, and generation of hydroxyl group. Moreover, the observation of a small number of CO2 in the oxidation products within the three types of UHMWPE models seemed to be consistent with the reported by-products of the irradiation PE in air 33 and during oxidation. 34 Correlating the variation of –OH groups (Figure 8) with that of small molecular chains (Figure 9a, b and c) threw light to anticipate that: (i) the cross-linking reaction of the oxidation products, which appeared as small molecular fragments, was attributed to the mutual collis ion between radicals; and (ii) there were basically two stages including oxidative decomposition and cross-linked reaction of UHMWPE under high-temperature and oxygen-enriched environment (see: Sec. 3.2).

Figure 9. Variation of the number of molecular fragments during the RMD processes: (a) RMD simulation result of UHMWPE-20; (b) RMD simulation result of UHMWPE-40; (c) RMD simulation result of UHMWPE-80.

3.4. The mechanism of oxidative reaction. Figure 10 shows the change in the number of functional groups to be 13


Page 21 of 28

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


identified in the processes of RMD simulation. Basically, it all (Figure 10a, b and c) have illustrated that: (i) functional groups like carbonyl, hydroperoxide, ketone, carboxylic acid, and aldehyde were all formed within the oxidative degradation stage when simulation time was shorter than 270 ps; and (ii) the oxidative groups like carbonyls and carboxylic acids, which were mainly observed in the stage of chemical cross-linking, were identified post of completing all the O2 reactions. The distribution of groups in the RMD simulation with the three types of UHMWPE models clearly demonstrated the existence of a common reaction order: forming of hydroperoxides (ROOH) in the 22

early stage of oxidation which was followed by the formation of carbonyl groups. A few of aldehydes were detected in the early stage of oxidative degradation although such aldehydes were completely absent in the final oxidation products. The fact that the rapid increase in ketone was always accompanied with the commencement of number decreasing in hydroperoxide suggested that the forming of ketone was mainly attributed to the decomposition of hydroperoxide.

Figure 10.Variation of the number of groups during the RMD processes: (a) RMD simulation result of UHMWPE-20; (b) RMD simulation result of UHMWPE-40; (c) RMD simulation result of 14



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 28

UHMWPE-80.

Figure 11 shows the progressive evolution of oxidative reaction mechanism of UHMWPE-20 in 3-D views. Generally, the mechanism of thermal-oxidation of the three types of UHMWPE models, as analyzed by RMD simulations, seemed agreeable with those relevantly experimental analyses available in literatures.

22,23,35-37

0.2 ps

80 ps

92 ps

94ps

120 ps

240 ps

600 ps

1000 ps

Figure 11. 3-D views of oxidative reaction mechanism of UHMWPE-20 from RMD simulation.

The results derived from the analysis of architecture, oxidation products and functional groups illustrated that the oxidation mechanism of UHMWPE could be as presented in Scheme of reactions (Figure 12). According to the mechanism of thermal oxidation process, the formation of alkyl macroradicals (R•) and hydroperoxides (as seen in Scheme, reactions 1 and 3 – Figure 12) were simultaneous ly taking place. The continuous formation of ketones, alcohols and H2 O were mainly attributing to (i) the decomposition of ROOH and (ii) the reaction of alcoxyl radical (RO•) or hydroxyl radical (•OH) to release an ion of H from the UHMWPE chains (see: Figure 12, reactions 5, 6, 7 and 8). Furthermore, the continuous formation of new secondary alkyl macroradicals kept the oxidation process on. The proceeding of mutually transformation between the carboxylic acids, alcohols and peroxy radicals (ROO•) was continuously on as long as sufficient O2 molecules were available within the reaction system.

15


Page 23 of 28

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


Figure 12. Scheme of reactions

4. CONCLUSIONS A novel computational methodology for modeling RMD s imulation of a representative polymer– UHMWPE by using a reactive force field of ReaxFF was described. Effect of molecular weight of atomic model on simulation results was also systematically analyzed through the discussion of diffusion coefficient, molecular architecture, thermodynamic changes, distribution of resulted products and the corresponding reaction mechanism. Data from the RMD simulations of UHMWPE-40 and UHMWPE-80 are correlated with available experimental results almost within acceptable level of good precision in terms of diffusion coefficient, resulted oxidation products and the corresponding oxidation mechanism. It thus confirms that the proposed structural models are considered to be representative to the reactive interactions within the UHMWPE systems, and can be used for understanding the atomic level of inter-molecular activities which is normally unobtainable by any means of simulations. It also sheds light for utilizing these computational approaches to investigate the process of micro-reaction for other macromolecular materials, which is surely one of our future research goals. The RMD simulations in this study allow several conclusions to be drawn. They are typically: (i) Comparing the values of diffusion coefficient D for the models of UHMWPE-40 D40 and UHMWPE-80 D80 with experimental data, it suggests that the two UHMWPEs giving exact capacity to allow O2 molecules penetrating into interior of the materials. But, the value of D20 deviated largely from the other calculated values of D40 and D80 , and 16



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 28

from the experimental data. Cross-linked reaction was usually considered as the second stage of thermal-oxidative reaction process. In re-polymerization reaction, observation of a higher cross-linking suggested the need to carry out RMD simulation for atomic model with larger molecular weight like UHMWPE-40 or UHMWPE-80. Successful RMD simulation for polymers required the atomic model with larger molecular weight, which has tendency to give (a) the exact diffusion coefficient with small molecules in polymer and (b) the significantly cross-linked reaction as reported in experimental works available in literatures. (ii) Although RMD simulations in this study consisted of UHMWPE-20, UHMWPE-40 or UHMWPE-80, the derived data of thermodynamics, products distribution and reaction mechanism seemed comparable to the experimental results. Simulated results suggested that the effect of molecular weight of atomic model on the precision level in RMD simulation of homopolymers like UHMWPE for thermodynamics, reaction mechanism and distribution of oxidation products was negligible. (iii) Selection of molecular weight for atomic model of macromolecules depends on the investigative category in RMD simulations. Generally, atomic model with larger molecular weight is preferred for predicting (a) the diffusivity of small molecules in polymer and (b) the changes of molecular architecture. However, results of RMD simulation for thermodynamics, products distribution, and reaction mechanism seemed to suggest smaller molecular weight of atomic model being sufficient to ensure the accuracy of simulation. Results in this study provide basic data and theoretical support for further RMD simulation to investigate chemical reaction of the interested macromolecules. AUTHOR INFORMATION Corresponding author *E-mail: [email protected] Notes The authors declare no competing financial interest. ACKNOWLEDGMENTS We are grateful for the financial support of the National Natural Science Foundation of China (No. 51275361) and the National Scientific and Technical Project of China (No. 2013CB632303). REFERENCES (1) Mueller, J.E.; van Duin, A.C.T.; Goddard, W.A. Application of the ReaxFF Reactive Force Field to Reactive Dynamics of Hydrocarbon Chemisorption and Decomposition. J. Phys. Chem. C 2010, 114(12), 5675-5685. (2) Chen, B.; Wei, X. Y.; Yang, Z. S.; Liu, C.; Fan, X.; Qing, Y.; Zong, Z. M. ReaxFF Reactive Force Field for 17


Page 25 of 28

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


Molecular Dynamics Simulations of Lignite Depolymerization in Supercritical Methanol with Lignite-Related Model Compounds. Energy Fuels 2012, 26(2), 984-989. (3) Nielson, K. D.; van Duin, A. C.; Oxgaard, J.; Deng, W. Q.; Goddard, W. A. Development of the Reaxff Reactive Force Field for Describing Transition Metal Catalyzed Reactions, with Application to the Initial Stages of the Catalytic Formation of Carbon Nanotubes. J. Phys. Chem. A 2005, 109(3), 493-499. (4) Mueller, J. E.; van Duin, A. C.; Goddard III, W. A. Development and Validation of Reaxff Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. C 2010, 114(11), 4939-4949. (5) Leininger, J. P.; Minot, C.; Lorant, F. Two Theoretical Simulations of Hydrocarbons Thermal Cracking: Reactive Force Field and Density Functional Calculations. J. Mol. Struct. 2008, 852(1), 62-70. (6) Salmon, E.; van Duin, A. C.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. Early Maturation Processes in Coal. Part 2: Reactive Dynamics Simulations Using te Reaxff Reactive Force Field on Morwell Brown Coal Structures. Org. Geochem. 2009, 40(12), 1195-1209. (7) Salmon, E.; van Duin, A. C.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. Thermal Decomposition Process in Algaenan of Botryococcus Braunii Race L. Part 2: Molecular Dynamics Simulations Using the Reaxff Reactive Force Field. Org. Geochem. 2009, 40(3), 416-427. (8) Nyden, M. R.; Stoliarov, S. I.; Westmoreland, P. R.; Guo, Z.X.; Jee, C. Applications of Reactive Molecular Dynamics to the Study of the Thermal Decomposition of Polymers and Nanoscale Structures. Mater. Sci. Eng. A 2004, 365(1), 114-121. (9) Stoliarov, S. I.; Westmoreland, P. R.; Nyden, M. R.; Forney, G.P. A Reactive Molecular Dynamics Model of Thermal Decomposition in Polymers: I. Poly (Methyl Methacrylate). Polym. 2003, 44(3), 883-894. (10) Stoliarov, S. I.; Lyon, R. E.; Nyden, M. R. A Reactive Molecular Dynamics Model of Thermal Decomposition in Polymers. II. Polyisobutylene. Polym. 2004, 45(25), 8613-8621. (11) Smith, K. D.; Bruns, M.; Stoliarov, S. I.; Nyden, M. R.; Ezekoye, O. A.; Westmoreland, P. R. Assessing the Effect of Molecular Weight on the Kinetics of Backbone Scission Reactions in Polyethylene Using Reactive Molecular Dynamics. Polym. 2011, 52(14), 3104-3111. (12) Jee, C. S. Y.; Guo, Z. X.; Stoliarov, S. I.; Nyden, M. R. Experimental and Molecular Dynamics Studies of the Thermal Decomposition of a Polyisobutylene Binder. Acta Mater. 2006, 54(18), 4803-4813. (13) Kemper, T. W.; Sinnott, S. B. Hyperthermal Atomic Oxygen and Argon Modification of Polymer Surfaces Investigated by Molecular Dynamics Simulations. Plasma Process. Polym. 2012, 9(7), 690-700. 18



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 28

(14) Chenoweth, K; Van Duin, A. C. T.; Goddard III, W. A. Reaxff Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112(5), 1040-1053. (15) Official Webs ite of Materials Studio. http://accelrys.com/products/collaborative-science/biovia-materials-studio/ (accessed Nov 7, 2015). (16) Kirkpatrick, S.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220(4598), 671-680. (17) Sun, H. COMPASS: an Ab Initio Force-Field Optimized for Condensed-Phase Applications Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102(38), 7338-7364. (18) Kantola, A.; Villar, H. O.; Loew, G. H. Atom Based Parametrization for a Conformationally Dependent Hydrophobic Index. J. Comput. Chem. 1991, 12(6), 681-689. (19) Ewald, P. P. Evaluation of Optical and Electrostatic Lattice Potentials. Ann. Phys. 1921, 369(3), 253-287. (20) Becke, A. D. a Multicenter Numerical Integration Scheme for Polyatomic Molecules. J. Chem. Phys. 1988, 88(4), 2547-2553. (21) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31(3), 1695-1697. (22) Bracco, P.; Brach del Prever, E. M.; Cannas, M.; Luda, M. P.; Costa, L. Oxidation Behaviour in Prosthetic UHMWPE Components Sterilised with High Energy Radiation in A Low-Oxygen Environment. Polym. Degrad. Stab. 2006, 91(9), 2030-2038. (23) Costa, L.; Carpentieri, I.; Bracco, P. Post Electron-Beam Irradiation Oxidation of Orthopaedic UHMWPE. Polym. Degrad. Stab. 2008, 93(9), 1695-1703. (24) Jost, S.; Biswas, P.; Schüring, A.; Kärger, J.; Bopp, P. A.; Haberlandt, R.; Fritzsche, S. Structure and Self-Diffusion of Water Molecules in Chabazite: a Molecular Dynamics Study. J. Phys. Chem. C 2007, 111(40), 14707-14712. (25) Daly, B. M.; Yin, J. Subsurface Oxidation of Polyethylene. J. Biomed. Mater. Res. 1998, 42(4), 523-529. (26) Pauly, S. Permeability and Diffusion Data; Brandrup, J., Immergut, E.H., Grulke, E.A., Abe, A., Bloch, D.R., Eds.; John Wiley & Sons: New York, U. S., 2003. (27) Rocha, M.; Mansur, A.; Mansur, H. Characterization and Accelerated Ageing of UHMWPE Used in Orthopedic Prosthesis by Peroxide. Mater. 2009, 2(2), 562-576. (28) Enomoto, I.; Katsumura, Y.; Kudo, H.; Soeda, S. Graft Polymerization Us ing Radiation-Induced Peroxides and Application to Textile Dyeing. Radiat. Phys. Chem. 2011, 80(2), 169-174. 19


Page 27 of 28

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


(29) Widmer, M. R.; Heuberger, M.; Vörös, J.; Spencer, N. D. Influence of Polymer Surface Chemistry on Frictional Properties Under Protein-Lubrication Conditions: Implications for Hip-Implant Des ign. Tribol. Lett. 2001, 10(1-2), 111-116. (30) Chen, J. S.; Zhu, F. Y.; Pan, H. C.; Cao, J. Q.; Zhu, D. Z.; Xu, H. J.; Cai, Q.; Shen, J. G.; Chen, L. H.; He, Z. G. Surface Modification of Ion Implanted Ultra High Molecular Weight Polyethylene. Nucl. Instr. and Meth. in Phys. Res. B 2000, 169(1), 26-30. (31) Shi, W.; Dong, H.; Bell, T. Wear Performance of Ultra High Molecular Polyethylene Ion Implanted Weight. Surf. Eng. 2003, 19(4), 279-283. (32) Chen, W.; Duan, H. T.; Hua, M.; Gu, K. L.; Shang, H. F.; Li, J. Comparison of Oxidation Resistance of UHMWPE and POM in H2 O2 Solution from ReaxFF Reactive Molecular Dynamics Simulations. J. Phys. Chem. B 2014, 118(34), 10311-10318. (33) Ivanov, V. S. Radiation Chemistry of Polymers; CRC Press: Florida, U. S., 1992. (34) Barabas, K.; Iring, M.; Kelen, T. Study of the Thermal Oxidation of Polyole-Fins. V. Volatile Products in the Thermal Oxidation of Polyethylene. J. Polym. Sci.: Polym. Sympo. 1976, 57(1), 65-71. (35) Costa, L.; Luda, M. P.; Trossarelli, L.; Del Prever, E. B.; Crova, M.; Gallinaro, P. In Vivo UHMWPE Biodegradation of Retrieved Prosthesis. Biomaterials 1998, 19(15), 1371-1385. (36) Brunella, V.; Bracco, P.; Carpentieri, I.; Paganini, M. C.; Zanetti, M.; Costa, L. Lifetime of Alkyl Macroradicals in Irradiated Ultra-High Molecular Weight Polyethylene. Polym. Degrad. Stab. 2007, 92(8), 1498-1503. (37) Medhekar, V.; Thompson, R. W.; Wang, A.; McGimpsey, W. G. Modeling the Oxidative Degradation of Ultra-High-Molecular-Weight Polyethylene. J. Appl. Polym. Sci. 2003, 87(5), 814-826.

20



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 Graphic

21


Page 28 of 28

Effect of Molecular Weight on Reactive Molecular Dynamics

Recommend Documents