Quantitative Comparison of Atomistic Simulations with Experiment for

Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

Quantitative Comparison of Atomistic Simulations with Experiment for a Cross-Linked Epoxy: A Specific Volume−Cooling Rate Analysis Ketan S. Khare*,†,‡ and Frederick R. Phelan, Jr.*,‡ †

Department of Physics, Georgetown University, 37th and O Streets, N.W., Washington, D.C. 20057, United States Material Science and Engineering Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899-8542, United States

‡

S Supporting Information *

ABSTRACT: Cross-linked epoxy thermosets, like all glassforming viscoelastic materials, show both a temperature and rate dependence in their thermomechanical properties. However, accounting for rate effects on these properties using molecular dynamics (MD) simulations and making quantitative comparison with experimental measurements has proven to be a difficult task due to the extreme mismatch between experimental and computationally accessible cooling rates. For this reason, the effect of cooling rate on material properties in glass-forming systems (including epoxy networks) has been mostly ignored in computational studies, making quantitative comparison with experimental data nebulous. In this work, we investigate a strategy for modeling rate effects in an epoxy network based on an approach that uses theoretically informed simulation and analysis protocols in combination with material specific time−temperature superposition (TTSP) data obtained from experimental measurements. To illustrate and test the strategy, we build and study an atomistic model of a cross-linked epoxy network. Molecular dynamics simulations are used to model the specific volume as a function of temperature across the glass transition from the rubbery to the glassy state using a total of five computationally accessible cooling rates. From the trends thus identified, we pinpoint the temperatures at which the models show rubbery and glassy behavior and use this information to calculate the values of the glass transition temperature (Tg) for each of the different cooling rates. Comparison with experimental data obtained from the literature (for the identical epoxy network) shows that our computations successfully predict the trends in specific volume in the rubbery and the glassy regions within 0.5%. We then compare the Tg values obtained from the data analysis with those calculated using the TTSP data obtained from the literature. Excellent agreement is found, and the Tg values from the two different methods are within 1.5% for all cooling rates. While our MD simulations do not replicate the experimental cooling rates, the agreement between these two sets of Tg values quantitatively relates the computational and experimental data sets. This agreement indicates that atomistic simulations can reliably capture the molecular mechanisms underlying viscoelasticity in this cross-linked epoxy (even when cooling rates differ by orders of magnitude) and that volume−rate analysis in conjunction with TTSP is a reliable method to computationally characterize certain classes of glass-forming materials. We believe that the general paradigm and protocols developed in this work should in principle be extensible as an efficacious means to integrate theory, computation, and experiment for the analysis of amorphous macromolecular materials.

1. INTRODUCTION Integration of experiment, computation, and theory for the study of commercially important materials is essential for connecting fundamental research with industrial applications; hence, such integration has been highlighted as one of the main goals of the U.S. Materials Genome Initiative (MGI).1 Crosslinked epoxy networks (or epoxy resins) comprise one of the most important classes of thermosetting polymers used in industry.2 Wide use of cross-linked epoxy networks results from a combination of versatile properties including mechanical strength and toughness, corrosion resistance, and desirable thermal, adhesive, and electrical characteristics.2 Even as the conventional uses of cross-linked epoxy networks continue to © XXXX American Chemical Society

grow, there has been a surge of interest in the use of crosslinked epoxy networks as host matrices for fillers in both composites and nanocomposites to improve mechanical, thermal, or electrical properties relative to the neat polymer networks.3−6 For such applications, the low viscosity of the mixture of the reactants before the gel point is advantageous for composite processingfor instance, the impregnation of fibers with cross-linked epoxy during the preparation of fiberreinforced composites.7 This surge of interest has in turn Received: June 20, 2017 Revised: December 13, 2017

A

DOI: 10.1021/acs.macromol.7b01303 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

with experimental data from the literature. Our approach is driven by the pursuit of the integration of experiment, computation, and theory, as envisioned in the goals of the MGI as a means to make more reliable materials with greater time and cost efficiency.1 For experiment, we extracted data from selected works32−35 on a specific material (described in the next section) in the literature. The key factor for our selection was the availability of data sets that could be directly compared with our simulations. Specifically, while it is far more common for experiments to report glass transition behavior using dynamic mechanical thermal analysis (DMTA) rather than by using volumetric trends using dilatometry, volumetric trends are almost exclusively17 used in atomistic simulations due to the ease of calculations. Since we seek direct quantitative comparison between simulation and experiment, only experimental works that present temperature and rate dependence of vsp were of direct interest in this case and other comparable methods of obtaining the Tg of cross-linked epoxy networks were excluded For computation, we performed MD simulations on atomistically detailed models of cross-linked epoxy. The molecular models were each created starting from a simulation of a molten mixture of the epoxy monomers and the cross-linkers. As will be explained in the subsequent section, we combined the simulated annealing (SA) method26,36,37 and the directed diffusion (DD) method18 to create fully cured and well-relaxed atomistic models of a cross-linked epoxy network. The volumetric behaviors of these models during cooling at five different cooling rates were then compared with experimental behavior obtained from the literature.32−35 For theory, we used the TTSP.29−31 Although the time scales accessible via atomistic simulations are short as compared to experiment, atomistic simulations can be performed at extremely high values of temperature.38 At sufficiently high values of temperature, vsp−T trends are expected to be rate independent, which is characteristic of rubbery behavior. Thus, even at the high cooling rates typically used in atomistic simulations, it should be possible to study the transition from rubbery to glassy behavior if simulations are performed at sufficiently elevated temperatures. Further, the values of Tg thus obtained can be quantitatively related to experimental values obtained at lower cooling rates using TTSP via the trend in the time-shift factor (aT) with the temperature (T). For linear viscoelastic behavior, the aT−T trend is unique for a given material irrespective of the viscoelastic phenomena being studied, which allows the creation of master curves for viscoelastic properties.30 This uniqueness of the aT−T trend is evidence that the molecular mechanisms underlying all linear viscoelastic phenomena are equivalent.32 Finally, wherever possible, we present our simulation data alongside experimental data obtained from the literature in the same or adjacent figures for emphasizing direct comparisons. In summary, although many reports of atomistic simulations of cross-linked epoxy networks are available in the literature,17,18 direct quantitative comparison has been relatively rare. The focus of our simulation effort here is to establish protocols that allow us to make direct quantitative comparison between simulation and experiment. In the subsequent section, we describe the computational methods and analysis techniques used in this work. This section is followed by the results and discussion section that has been organized with subsections based on our key findings. Finally, we close our paper with a conclusion section that presents a summary and also a list of

raised fundamental research questions about the nonintuitive physics operating in the interphase region of such composites, which arises due to the matrix−filler interactions.8−14 The inclusion of specific chemical interactions in atomistically detailed molecular dynamics simulations (i.e., atomistic simulations) makes them well-suited for gaining insight at such molecular length scales.10,15−17 An extensive review of molecular scale simulations of thermosetting polymers was recently published.17 A large fraction of the review17 was focused on the simulation of atomistic models of cross-linked epoxy networks. This review documents all of the significant advances that have recently occurred in the field.17 However, as can be seen in that review17 and elsewhere,18 the extents of curing/conversion of the model networks in a majority of the simulation studies are far from completion.17,18 By contrast, the experiments in the literature19−25 show that (1) the extent of conversion of the network can be driven essentially to completion (conversion was measured using both infrared spectroscopy 21−24 and differential scanning calorimetry19,20,24,25) and that (2) the glass transition temperature (Tg) and the viscoelastic behavior of the network are sensitive19−25 to the extent of conversion. Partly because of this sensitivity to the extent of conversion, quantitative agreement between such simulations and experiments has been a challenge to document.14 R. Khare and co-workers18,26−28 have performed simulations of fully cured atomistic models for cross-linked epoxy networks of many different chemistries. In their works,18,26−28 it has been consistently shown that (1) the network densities below Tg, which were obtained via simulation are about 2% lower than those obtained via experiment for the same network, and (2) the values of Tg obtained via simulation are about 30−110 K higher than those obtained from experiment. Both of these differences were expected due to the high cooling rates used in the simulations compared to those typically used in experiments. Further, the experimental and computational values of Tg showed good agreement when the differences in the cooling rates were taken into account using the time−temperature superposition principle (TTSP).29−31 However, in those works,18,26−28 only one cooling rate was used to relax the network structure of the model. Further, the trend in the specific volume (vsp) with temperature (T) was not directly compared with experimental data. In atomistic simulations of thermosets, the characteristic vsp−T trend during cooling is the phenomena most frequently used to observe the glass transition and calculate the value of the Tg.17 The dependence of both the Tg and the vsp−T trend on the cooling rate arises because of the underlying molecular mechanisms for viscoelasticity. The selection of a particular cooling rate effectively creates an observation window for the study of the relaxation of vsp as a function of T. The molecular mechanisms underlying the temperature and cooling rate trends of specific volume are the same as those governing the response of the material to mechanical deformations, such as shear, elongation, or creep.32 Thus, the ability of simulations to successfully capture (qualitatively and more importantly quantitatively) both the temperature and the rate dependence of vsp of a thermoset is an indicator that the mechanisms governing predictions of Tg, elastic moduli, shear moduli, creep compliances, etc., are also properly modeled. In this work, we use atomistic simulations of cross-linked epoxy to make direct quantitative comparisons of temperature and rate dependent trends of vsp and other related quantities B


Article

Macromolecules recommendations for continuing progress toward the integration of experiments, computation, and theory for glass-forming polymeric networks.1

Table 1. Simulation Details and Parameters

2. COMPUTATIONAL METHODOLOGY 2.1. Network Chemistry. An amine-cured epoxy resin network was studied here. The industrial production of the class of epoxy monomers based on the diglycidyl ether of bisphenol A (DGEBA) chemistry results in a mixture of oligomers with varying degrees of polymerization (n). Here, we use DGEBA with n = 2 to represent the commercial product Epon 1001F.39 4,4′-Diaminodiphenyl sulfone (4,4′-DDS) was used as the cross-linker. The chemical structures of both molecules are shown in Figure 1.40 The volumetric and glass transition properties of this epoxy network were the subject of extensive experimental investigations in previous work.32−35

detail/parameter

description/value

force field partial charges short-range cutoff long-range van der Waals long-range Coulombic constraint for bonds and angles containing hydrogen simulation time step temperature (T) pressure (P) Nosé−Hoover Tdamp57−59 Nosé−Hoover Pdamp57−59 box size

GAFF45−47 AM1-BCC48,49 9Å tail pppm54 RATTLE55,56 1 fs varies 5 MPa 100 fs 1000 fs 13 nm

2.3. Structure Preparation. Five replicas of the crosslinked epoxy network were created from stoichiometric mixtures of Epon 1001F (monomer) and 4,4′-DDS (crosslinker) molecules as follows. Two molecules of the monomer (Figure S1 of the Supporting Information)40 and one molecule of the cross-linker were added to a simulation box and replicated along the Cartesian axes such that the resulting simulation box was roughly cubic and contained 1458 monomer molecules and 729 cross-linker molecules in a lattice arrangement. Thereafter, this box was simulated to form molten mixtures at a temperature of 800 K and pressure of 10.1 MPa (100 atm) for a simulation time of about 10 ns. We note that the chosen simulation temperature is significantly above the experimental melting points33 of Epon 1001F and 4,4′-DDS (about 350 and 450 K, respectively), and this temperature condition facilitated the homogeneous mixing of the molecules in the simulation box within time scales of a few nanoseconds. The five replicas were generated using different seeds for the random number generator of the initial velocities. The distribution of the molecules in each of the five replicas diverged, and thus, each of the five replicas is effectively independent.38 The curing/polymerization of the mixture was performed using the Simulated Annealing37 (SA) algorithm. SA was initially developed as an effective algorithm to provide reasonable solutions to the traveling salesman problem (TSP).37 The algorithm was previously adapted to create atomistically detailed models of polystyrene36 and more recently to create models of cross-linked epoxy networks.26,28 Here, we used custom code to implement the SA method, which provided us with a connectivity sequence between the monomers and the cross-linkers. Using this connectivity sequence, the monomers and cross-linkers are allowed to diffuse using the Directed Diffusion18 (DD) strategy. The purpose of DD is to facilitate the diffusion of reacting partners of the monomers and the cross-linkers using MD simulations at an accelerated time scale of a few nanoseconds rather than the hours that are used to cure epoxy networks in experiments.34 Further information about the SA and DD methods is available in the relevant literature18,26 and also in the Supporting Information. The appropriate angles, dihedrals, and impropers were added to account for the newly created bonds according to GAFF.45−47 Partial charges were also updated using AM1BCC.48,49 Finally, the five replicas of the cross-linked epoxy model structures were allowed to relax using MD simulations at a temperature of 820 K for a simulation time of 10 ns.

Figure 1. Chemical structure of (a) Epon 1001F and (b) 4,4′-DDS.

2.2. Simulation Details. Molecular dynamics (MD) simulations of atomistically detailed models were used in this work. The details of the methods used here are almost identical to those in previous work10,15,18,26−28,41−44 but have been included here for the sake for completeness. The all-atom general Amber force field45−47 (GAFF) was used to describe the molecular models. Partial charges were calculated using the Austin Model 1 with bond charge correction model (AM1BCC).48−50 The AmberTools 16 software package was used to assign the force field and the partial charges.51 Additional details are available in the Supporting Information. All simulations were performed using the LAMMPS52 simulation package. The pairwise van der Waals interactions were truncated at a distance of 9 Å, and a tail correction53 was applied to the energy and pressure of the system to account for the residual. The pairwise Coulombic interactions were calculated directly up to a distance of 9 Å. Long-range Coulombic interactions were computed in reciprocal space using the particle-particle particle-mesh (pppm) method.54 The RATTLE55,56 algorithm was used to constrain the bonds and angles containing hydrogen atoms. Simulations were performed using a time step of 1 fs. We used the Nosé−Hoover57−59 thermostat and barostat to maintain isothermal−isobaric (constant number of particles, pressure, and temperature or constant NPT) conditions. The fluctuation times of the thermostat and barostat (Tdamp and Pdamp) were 100 and 1000 time steps, respectively. Unless otherwise specified, the value of the pressure used for the simulations was 5 MPa (49.35 atm), corresponding to that used to obtain experimental data in the literature.32,34 We note that the pressure used in our work was not varied in order to reach any target values of density or vsp. The value of vsp at a given temperature and cooling rate is solely the result of the applied force field, simulation protocol, and the network topology of the simulated model. These simulation details have also been summarized in Table 1. C


Article

Macromolecules

Table 2. Details of Various Cooling Rates Used in Simulations and Experimental Values Extracted from Literature Dataa Tg (K)

temperature range (K)

a

name

cooling rate (K s−1)

cooling rate (convenient units)

slow−1× medium−3× fast−9× very fast−30× ultrafast−100× Plazek and Choy35 (PC) Bero and Plazek32 (BP1) Bero and Plazek32 (BP2) Bero and Plazek32 (BP3) Bero and Plazek32 (BP4) WLF reference35

5.556 × 109 1.667 × 1010 5.000 × 1010 1.667 × 1011 5.556 × 1011 8.333 × 10−2 5 × 10−5 8.33 × 10−4 4.17 × 10−3 1.5 × 10−2 2.64 × 10−3

5.556 GK s−1 16.67 GK s−1 50.00 GK s−1 167.7 GK s−1 555.6 GK s−1 5 K min−1 0.003 K min−1 0.050 K min−1 0.250 K min−1 0.90 K min−1 0.157 K min−1

glassy

rubbery

280 to 380 620 to 720 280 to 380 660 to 760 280 to 380 710 to 810 280 to 380 730 to 830 280 to 380 745 to 845 not applicable; experimental data extracted from literature; uncertainty estimates not available

calculated

WLF

489 ± 1.0 501 ± 1.8 516 ± 3.1 536 ± 0.6 556 ± 1.5 407.3b 399.3 402.1 403.7 404.9 403.2

492 502 514 529 548 407 399 402 404 405

Uncertainty is not available for data extracted from the literature. bCalculated from extracted trend in specific volume with temperature.

2.4. Cooling Model Structures across the Glass Transition. Each of the five replicas of the cross-linked epoxy models was cooled from a high temperature of 820 K in the rubbery state to a low temperature of 140 K in the glassy state using temperature steps of 5 K. Five different cooling rates were used for each replica, which are summarized in Table 2. The different cooling rates were achieved by varying the simulation time at each temperature step from as long as 900 ps (slow−1×) to as short as 9 ps (ultrafast−100×). Data sets for analysis were collected during the last half of each temperature step. We note that even the slowest cooling rate used here is more than 10 orders of magnitude faster than those typically used in experiments.32−35 We also note that the slow−1× cooling rate used in this work is similar, albeit slightly slower, than that used in previous work.18,26,28 In Figure 2a, representative trends of specific volume with temperature at two cooling rates (fast and slow also denoted by superscript 1 and 2, respectively) are shown. Using the trend for the slow cooling rate as a reference, the relative specific volume (Δvsp) can be calculated as a difference between the two trends. Δvsp = vsp − vspref

(1)

The relative specific volume calculated in this way, which is shown in Figure 2b, allows an objective identification of T1, which is the temperature at which the specific volume has insufficient time to relax to equilibrium and hence deviates from rubbery behavior. Above T1, the trends are cooling rate independent, and hence Δvsp is zero. Thus, a linear fit to the vsp trend above T1 can be used to obtain the equilibrium rubbery line. In Figure 2b, it can be seen that below a temperature of T1 the relative specific volume (Δvsp) increases from zero and reaches a constant value at temperature T2. In the literature, this temperature T2 has been shown to be cooling rate independent, in contrast to the rate dependent temperature T1.32 Using a linear fit to the trend below T2, the glassy line can be obtained. While the discussion of the representative trends in Figure 2 involves two cooling rates, we used five different cooling rates for the simulations in our work. 2.5. Calculation of Tg. We used the vsp−T trends for each of the five replicas of the model structures at each of the five cooling rates to calculate the glass transition temperature (Tg). Specifically, the rubbery and the glassy regions for each of the five cooling rates were identified using the Δvsp−T trend, and

Figure 2. Schematic of the temperature trend of specific volume and two derived quantities at two different cooling rates. (a) Specific volume (vsp). (b) Relative specific volume using the slow cooling trend as a reference (Δvsp). (c) Departure of specific volume from the equilibrium rubbery line (above Tg) and glassy line (below Tg). As the temperature is reduced, the trend in vsp deviates from the equilibrium rubbery line at temperature T1. After the glass transition at temperature Tg, the trend merges with the glassy line at T2. The superscripts refer to the two cooling rates: 1 is the fast cooling rate; 2 is the slow cooling rate.

D


Article

Macromolecules

Figure 3. Specific volume (vsp) versus temperature (T). The simulation trends show cooling rate independence at high temperatures and gradually diverge from the rubbery trend in a cooling rate dependent manner. The extrapolated rubbery trend shows excellent agreement with the experimental trend extracted from the literature.35 Error bars have been removed for readabilit but are shown in Figures 4 and 5. Trends in regions A−C are shown in Figures 4−6, respectively.

the Tg was designated to be the temperature at which the linear fits to the vsp−T trend in the rubbery and glassy states intersected. A summary of the temperature ranges used to calculate Tg is presented in Table 2. In the case of the two fastest cooling rates (i.e., very fast−30× and ultrafast−100×) where the rubbery state extended above the highest simulated temperature, we extrapolated the rubbery trend obtained from the slowest cooling rate to higher values of temperature. Using the glass and the equilibrium rubbery line, as references, we can also calculate the departure of the specific volume δ(vsp) using eq 2: δ(vsp) = vsp − v line

cally equivalent. We prefer to use the WLF relation (eq 3) in our work: log aT = log

qr qT

=

−C1(T − Τr) C2 + (T − Tr )

(3)

where the WLF parameters are C1 = 19.26, C2 = 50, and Tr = 403.2 K.35 Further, qr and qT are the values of cooling rates that would cause the network to undergo glass transition at temperatures Tr and T, respectively. Using literature32 data for the dependence of Tg on cooling rate for the same system (Table 2), we estimate that the value of the cooling rate (qr) at the reference temperature T = Tr is 2.64 × 10−3 K s−1. We note that the trend in the time-shift factors (aT) with temperature for the same material using volume−rate dependent processes is also available in the literature.32 2.7. Data Extraction and Uncertainty Quantification. In order to make quantitative comparisons with experiments in the literature, data sets were extracted from various works in the literature.32−35 For this purpose, the WebPlotDigitizer61 software was used to extract data from figures. Unless otherwise stated, “error bars” in figures or uncertainty associated with quantities in the text reflect the “standard error of mean” for the five model structures. Where error bars are not shown, unless otherwise stated, the size of the error bars is similar to or smaller than the size of the symbols.

(2)

where vline is the value of the equilibrium rubbery or glassy line above or below Tg, respectively. The quantity δ(vsp) has been used in previous experimental work to study the breadth of the glass transformation region.32 Representative δ(vsp)−T trends are shown in Figure 2c. As discussed earlier, it can be seen that while the values of T1 are cooling rate dependent, the values of T2 are cooling rate independent. Thus, the use of the trends in Δvsp and δ(vsp) with temperature allow us to identify the rubbery and the glassy regions with a greater degree of objectivity than that possible using the vsp−T trend alone. In the literature, it was been shown that significant uncertainty can be introduced by a subjective identification of the rubbery and the glassy regions.60 2.6. WLF and TTSP. The Vogel−Tammann−Hesse− Fulcher30 (VTHF or VFT) parameters for cross-linked epoxy network formed by Epon 1001F and 4,4′-DDS are available in the literature.35 The WLF and VTHF equations are mathemati-

3. RESULTS AND DISCUSSION 3.1. General Features of the Volumetric Trends. The average trends in the specific volume (vsp) of the five replicas at five different cooling rates with temperature (T) are shown in E


Article

Macromolecules

temperature; the trends appear to be vertically shifted with respect to each other (region B in Figure 3; cf. Figure 5). (5) The experimental trend of vsp in the vicinity of the Tg shows a visible kink; however, such an abrupt transition is not seen for the simulation trends. (6) Extrapolation of the high temperature rubbery trend obtained from the simulations (to lower temperature) converges with the experimental trend in the rubbery state. 3.2. Specific Volume Trend in the Rubbery State Is Independent of Cooling Rate. The vsp−T trends at high temperature are shown in Figure 4 (cf. region A of Figure 3). As can be seen from the figure, the specific volume at a given temperature is independent of the cooling rate with the partial exception of the trend for the ultrafast−100× cooled models. This behavior is characteristic of the rubbery state where the network has sufficient time to relax to the corresponding thermodynamic state point at each step. The trend in the ultrafast−100× cooled models shows a slight deviation from the equilibrium rubbery line at temperatures lower than 780 K, indicating that below this point the short simulation time (9 ps) at each temperature step is insufficient to obtain the relaxation of vsp to the state point value. 3.3. Specific Volume Trends in the Glassy State Are Cooling Rate Dependent. In contrast to the high temperature observations, the behavior of the specific volume at low temperature shows a distinct cooling rate dependence. In Figure 5, the low temperature trends in specific volume are shown (cf. region B of Figure 3). This shows that the value of the specific volume at a given temperature increases with an increase in the cooling rate. This effect is even more pronounced for the literature data in which case the experimental vsp values are about 2% lower than those obtained via simulation. 3.4. Specific Volume Can Be Accurately Predicted Using Simulations. Since the rubbery state values of the specific volume and temperature at a given pressure represent thermodynamic state points that are independent of cooling rate, it is possible to make quantitative comparisons of the simulation results with the experimental data. However, the experimental data set is only available up to a temperature of about 460 K,35 which is significantly below the values of the glass transition temperature at the cooling rates used in the simulations. Thus, to make a state-point comparison, we use an extrapolation of the high temperature equilibrium simulation data to the lower temperature range spanned by the experimental data. This is enabled and justified by the aforementioned cooling rate independence of the specific volume in the rubbery state. By combining this extrapolation with a shift of the low temperature data, we show the entire experimental glass transition event can be matched to great accuracy. First, to construct the extrapolation, we assume that the value of the coefficient of volumetric thermal expansion (αv) is constant in the equilibrium state. For the equilibrium state, we use the rubbery state simulation data in the temperature range 700−820 K, at the slow−1× cooling rate. The thermal expansion coefficient is defined by the relationship

Figure 3. Details of the different cooling rates are summarized in Table 2. Also included in the figure is the vsp−T trend obtained from the experimental work of Plazek and Choy.35 The following general observations can be made from the figure, which will discussed in greater detail in subsequent subsections: (1) At high temperatures, all the trends obtained via simulation for the five different cooling rates converge onto the same high temperature rubbery trend (region A in Figure 3; cf. Figure 4). (2) As the system is cooled, but still above Tg, the

Figure 4. Specific volume (vsp) versus temperature (T) in the rubbery state. This figure is region A in Figure 3. With the partial exception for the fastest rate, the trends show cooling rate independence.

rubbery region trends begin to diverge starting with the trend obtained at the fastest cooling rate and in order down through the slowest rate (region B in Figure 3; cf. Figure 5). (3) The experimental trend of vsp obtained from the literature data lies significantly below all the trends obtained via simulation (region C in Figure 3; cf. Figure 6a). (4) At temperatures below Tg, all the observed trends (whether from simulation or experiment) appear to have identical slopes at a given

Figure 5. Specific Volume (vsp) versus temperature (T) in the glassy state. This figure corresponds to region B in Figure 3. Unlike the rubbery state, the trends show definitive cooling rate dependence. It can be observed that all the slopes are nearly identical. The differences in the cooling rate span 12 orders of magnitude.

αv = F

1 dvsp vsp dT

(4) DOI: 10.1021/acs.macromol.7b01303 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 6. (a) Specific volume (vsp) versus temperature (T) showing the glass transition. This figure is region C in Figure 3. Error bars have been removed for readability but are shown in Figures 4 and 5. (b) Comparable experimental data extracted from Bero and Plazek.32 The dashed gray line shows the equilibrium rubbery line used by the authors to calculate the Tg.

Using the equilibrium state data along with eq 4, we find that the rubbery equation of state can be described by the relationship ln(vsp) = 0.000541T − 0.385

(5)

This equation of state is plotted in Figure 6a as the black line with long dashes (extrapolated rubbery trend). It can be seen that the extrapolation using eq 5 shows excellent agreement (within 0.3%) with the experimental data obtained from the literature.35 We can further attempt to match the experimental data using the empirical observation noted in Figure 5 namely, that the glassy trend in vsp with T shows a constant slope independent of cooling rate (for both simulation and experiment32) below a temperature of 390 K. If we attempt to match the experimental Tg value of 407.3 K by shifting the simulation trends vertically downward until they both overlap and intersect with the extrapolated trend, i.e., eq 5, we find that the predicted values of vsp are within 0.5% of the experimental values as shown in Figure 6a. The shifted data are reflected in the curve labeled “predicted glassy trend”. In totality, we see that the combination the high temperature extrapolation with this low temperature shift allows us to virtually match the entire experimental glass transition curve. We are not aware of a construction of this kind having been previously reported in the literature. The simulation results may also be qualitatively compared with experimental cooling rate data. In Figure 6b, we show the cooling rate dependence of the experimental vsp−T trend obtained from the literature.32 Comparison of Figures 6a and 6b shows a striking qualitative resemblance of the simulation and experimental data. This strongly indicates that the modeling captures the same underlying viscoelastic processes that we know are reflected in the experimental trend, despite the vast differences in cooling rates. 3.5. Relative Specific Volume Trends Aid Identification of the Rubbery and the Glassy State. The trend in the relative specific volume (Δvsp) vs T for the five different cooling rates is shown in Figure 7. These values were calculated using eq 1 with the slow−1× simulation cooling rate data as the reference trend. Additionally, we plot values obtained using the literature data35 relative to the same reference trend. All the

Figure 7. Relative specific volume (Δvsp) using the slow−1× cooling rate as a reference. Representative error bars (see section 2.7) have been shown. The values of Δvsp at Tg are as a guide. The simulation and experimental data both plateau at approximately the same value of T2 .

curves show similar features. At temperatures above the T1 corresponding to the cooling rate, Δvsp is zero (indicating departure from the equilibrium rubbery line). Below T1, the values of Δvsp begin to grow and become constant below a temperature of T2. The simulation data show a positive deviation relative to the reference point, while the experimental data have a negative deviation. We find the value of T2 obtained via the simulation data is 390 K and is independent of cooling rate. This is in excellent agreement with the experimental result (T2 = 391 K) reported in the literature.32 3.6. Value of Tg Increases with an Increase in Cooling Rate. Qualitatively, it can be inferred from the data already presented that the values of the Tg decrease with a decrease in cooling rates. Quantitatively, Tg values were calculated from the intersection of the linear fits to the temperature dependence of vsp in the rubbery and glassy regions, respectively. The Tg values for each cooling rate are reported in Table 2 along with the ranges for the rubbery and glassy states. In these calculations, G


Article

Macromolecules

Figure 8. (a) Departure of specific volume (δ(vsp)) versus temperature (T) in simulations compared to one experimental trend.35 Error bars (see section 2.7) for the simulation data are the size of the symbols. (b) Comparable experimental data extracted from Bero and Plazek.32

eq 5 was used to represent the equilibrium rubbery line for the very fast−30× and ultrafast−100× rates. As can be seen in Table 2, the value of Tg consistently increases with an increase in the cooling rate. The computational Tg values are on the order of 90−150 K higher than the experimental values obtained at much lower cooling rates. As we will discuss later, we find that the dependence agrees completely with expectations based on TTSP. 3.7. Breadth of the Glass Transformation Region Increases Consistently with an Increase in Cooling Rate. The simulation data for the departure of the specific volume δ(vsp), calculated using eq 2, are shown in Figure 8a. Comparable experimental data sets that were extracted from Bero and Plazek32 are shown in Figure 8b. The former figure also shows the experimental trend from Plazek and Choy35 for the purpose of comparison. As discussed earlier and shown in Figure 2c, the quantity δ(vsp) visualizes the breadth of the glass transformation region and also helps identify the rubbery and the glassy states. The figure shows that for both the computational and experimental data the breadth of the glass transformation region increases with an increase in cooling rate. Further, both sets of data clearly illustrate the cooling rate dependence of T1 and the cooling rate independence of T2. What is striking is that the computational glass transformation regions extend for well over 100 K in all cases, while the experimental regions are far narrower. This largely explains the lack of a characteristic kink in the vsp−T trends obtained from simulation as compared to experiment where it is nearly always evident. 3.8. Values of the Coefficient of Volumetric Thermal Expansion at Tg Are Cooling Rate Independent. The values of the coefficient of volumetric thermal expansion (αv) vs T as a function of cooling rate are shown in Figure 9. These values were calculated using a quartic polynomial fit to the temperature trends for vsp and eq 4. As can be seen in the figure, the glass transformation range for the simulation is much larger than that for the experiment. The peak αv values for both data sets are of the same order of magnitude, but the simulation values show a clear cooling rate dependence. Below Tg, this dependence begins to lessen, and then both the experimental and simulation data coincide below 390 K, which is the approximate value of T2 for both cases. We note that the simulation values show a gradual decrease with temperature as compared to the experimental data which shows a precipitous

Figure 9. Coefficient of volumetric thermal expansion (αv) versus temperature (T). The simulation trends show strong cooling rate dependence. But, the values of αv at their respective values of Tg for both simulation and experiment are approximately constant, despite drastic differences in the individual values of cooling rate and Tg. Representative error bars (see section 2.7) are shown.

drop (as might be reflective of a second-order phase transition). Interestingly, all the simulation values and the experimental value of αv are approximately the same value at their respective glass transition temperatures, despite the vast cooling rate and individual Tg differences. 3.9. TTSP Allows the Bridging of Material Behavior across More Than 10 Orders of Magnitude in Time Scale and 150 K in Temperature. In Figure 10, we plot the timeshift factors aT with temperature from three distinct sources: (1) experiments by Bero and Plazek,32 (2) the WLF equation (eq 3) using material specific parameters from the literature,35 and (3) the present work using the Tg values obtained from the cooling rate trends in simulations. As can be seen in the figure, the aT values obtained using the simulation data fall virtually on the same line in excellent agreement with the WLF trend from the literature35 and quantitatively relates the two data sets. This agreement demonstrates that despite the vast differences in cooling rate between experiment and simulation, the underlying viscoelastic processes are equivalent in both cases and are reliably captured by the atomistic simulation. Further, this H


Article

Macromolecules

the utility of TTSP and the WLF relationship for connecting viscoelastic behavior across vast differences in values of temperature and cooling rates. (6) The general Amber force field (GAFF)45,46 used here was originally parametrized for small organic molecules, not polymers, and we did not reparametrize it for our material. Nonetheless, our results show that it is able quantitatively predict both the rate and temperature dependent trend in vsp for this epoxy system. In addition to atomistic simulations of cross-linked epoxy,10,15,18,26−28,41−44 GAFF has been previously used to simulate lipid bilayers,62 ionic liquids,63 and hydration of organic molecules.64 Thus, GAFF appears to show a good degree of transferability. (7) As described earlier,18 our work demonstrates the importance of obtaining complete conversion for the reaction between monomers and cross-linkers for leveraging the predictive abilities of atomistic simulations. (8) Molecular models of the size used in this work (about 200 000 atoms) appear to capture glassy and rubbery state behavior as well as the glass transition. On the basis of our work, we have the following recommendation for researchers seeking quantitative agreement between experiments and simulations. (1) The use of multiple cooling rates and an analysis of their relative volumetric trend and departure of the specific volume are helpful in pinpointing the onset and end of the glass transformation. The equilibrium rubbery and glassy lines obtained in this manner provide a more robust method for calculating the value of Tg. (2) The availability of experimental shift factors at higher temperatures would aid the further integration of experiments and computation. As seen in Figure 10, there is a gap of about 80 K in temperature and 6 orders of magnitude in time scale in the trend between shift factor and temperature. We also note the following caveats to the methods that may affect the extension of the protocols in this work to future work for other properties or materials: (1) We specifically selected the matrix chemistry based on the availability of extensive experimental data sets on (a) cooling rate dependence of specific volume, (b) time-shift factors, and (c) the materialspecific WLF parameters. Such data sets (especially (a)) are not always available for various materials. The absence of the value of the cooling rate independent glass temperature T2 obtained using experiments is likely to introduce some additional subjectivity in determining the temperature range for the glassy lines. As a result, the degree of agreement possible between simulations and experiments may suffer due to the absence of such experimental data sets. (2) The models created here are free from defects or local heterogeneities. Such idealistic models have not affected the agreement between experiments and simulations. However, such differences between models and experimental samples can hamper integration between simulations and experiments of other properties or other networks. Further, the selective absorption of chemical species at interfaces in nanocomposites can perturb the network structure in the matrix−filler interphase region. Approaches such as the reactive MD (RMD) simulations may have to be considered to account for such effects in the simulation models.12−14 (3) An overlap in the α- and β-relaxation causes thermorheological complexity that precludes the straightforward use of the time−temperature superposition. Such an effect is especially true for material exhibiting strong β-relaxation.65 Further, the temperature window of β-relaxation effects also shifts to higher temperatures at faster rates typically used in

Figure 10. Shift factor (aT) versus temperature (T). The cooling rate (qT) equivalent to shift factors is shown on the right axis. We find excellent agreement between the WLF equation parametrized using experiments and the simulations. Error bars (see section 2.7) for the simulation Tg are the size of the symbols.

illustrates the utility of using TTSP in conjunction with volume−rate analysis to quantitatively relate these disparate data sets by bridging very large differences in time scales. All these methodssimulation, analysis based on theory, and TTSPused together provide a powerful means for integrating experiment, computation, and theory in the study of epoxies and other similar classes of glass-forming materials.

4. CONCLUSIONS In this work, we perform atomistic simulations of cross-linked epoxy and directly compare our results with those obtained using experiments in the literature. (1) The calculations predict that the vsp−T trend at high temperatures in the rubbery state is cooling rate independent, in agreement with experimental data for the same system. 32 Further, we show that using extrapolation of the equilibrium rubbery line obtained from the simulations to the rubbery region spanned by the experimental data, we can predict the experimental values of vsp to within 0.3% (Figure 6). The fact that this works despite the 250 K extrapolation indicates that the molecular mechanism underlying rubbery behavior is insensitive to both cooling rate and temperature and is well captured by the MD simulation. It also indicates that the assumption of a constant value of αv in the rubbery state is reasonable. (2) Using a reconstruction in which we shift the constant slope glassy lines to intersect with the extrapolated rubbery trend, we find that we can predict the value of vsp in the glassy state to within 0.5% (Figure 6). This agreement supports the observation made in the literature32 that the glassy lines are reached at the same temperature despite the drastic differences in the cooling rate. (3) For this epoxy system, the calculations predict that the glass transformation range ends at a temperature of T2 = 390 K, which agrees with the experimental32 value of T2 = 391 K (Figures 7 and 8). (4) The simulation values and the experimental value of αv are approximately the same at their respective glass transition temperatures, despite the vast cooling rate and individual Tg differences (Figure 9). (5) The trend in Tg with cooling rate in simulations shows excellent quantitative agreement with expectations from the WLF relationship (Figure 10 and Table 2). Such an agreement is a powerful demonstration of I


Article

Macromolecules simulations, as seen in previous work by one of us.28 (4) While GAFF45,46 is shown to work successfully for the models simulated in this work (and also works in the literature for cross-linked epoxy thermosets of a variety of chemistries18,26−28,44), the force field has not been designed or validated to operate for polymeric chemistry nor at the exceptionally high values of temperatures often used in atomistic simulations. In future work, we intend to perform creep deformation on these models to investigate the predictive abilities of atomistic models for mechanical properties.32

■

(7) Liu, Y.; Yang, J.-P.; Xiao, H.-M.; Qu, C.-B.; Feng, Q.-P.; Fu, S.-Y.; Shindo, Y. Role of Matrix Modification on Interlaminar Shear Strength of Glass Fibre/Epoxy Composites. Composites, Part B 2012, 43, 95− 98. (8) Putz, K. W.; Palmeri, M. J.; Cohn, R. B.; Andrews, R.; Brinson, L. C. Effect of Cross-Link Density on Interphase Creation in Polymer Nanocomposites. Macromolecules 2008, 41, 6752−6756. (9) Rittigstein, P.; Priestley, R. D.; Broadbelt, L. J.; Torkelson, J. M. Model Polymer Nanocomposites Provide an Understanding of Confinement Effects in Real Nanocomposites. Nat. Mater. 2007, 6, 278−282. (10) Khare, K. S.; Khabaz, F.; Khare, R. Effect of Carbon Nanotube Functionalization on Mechanical and Thermal Properties of CrossLinked Epoxy−Carbon Nanotube Nanocomposites: Role of Strengthening the Interfacial Interactions. ACS Appl. Mater. Interfaces 2014, 6, 6098−6110. (11) Grady, B. P. Effects of Carbon Nanotubes on Polymer Physics. In Carbon Nanotube−Polymer Composites; John Wiley & Sons, Inc.: 2011; pp 119−189. (12) Langeloth, M.; Sugii, T.; Böhm, M. C.; Müller-Plathe, F. Formation of the Interphase of a Cured Epoxy Resin Near a Metal Surface: Reactive Coarse-Grained Molecular Dynamics Simulations. Soft Mater. 2014, 12, S71−S79. (13) Langeloth, M.; Sugii, T.; Böhm, M. C.; Müller-Plathe, F. The Glass Transition in Cured Epoxy Thermosets: A Comparative Molecular Dynamics Study in Coarse-grained and Atomistic Resolution. J. Chem. Phys. 2015, 143, 243158. (14) Farah, K.; Leroy, F.; Müller-Plathe, F.; Böhm, M. C. Interphase Formation during Curing: Reactive Coarse Grained Molecular Dynamics Simulations. J. Phys. Chem. C 2011, 115, 16451−16460. (15) Khare, K. S.; Khare, R. Effect of Carbon Nanotube Dispersion on Glass Transition in Cross-Linked Epoxy−Carbon Nanotube Nanocomposites: Role of Interfacial Interactions. J. Phys. Chem. B 2013, 117, 7444−7454. (16) Khare, K. S. Thermo-mechanical Properties of Cross-linked Epoxy based Systems: A Molecular Simulation Study; Texas Tech University: 2013. (17) Li, C.; Strachan, A. Molecular Scale Simulations on Thermoset Polymers: A Review. J. Polym. Sci., Part B: Polym. Phys. 2015, 53, 103− 122. (18) Khare, K. S.; Khare, R. Directed Diffusion Approach for Preparing Atomistic Models of Crosslinked Epoxy for Use in Molecular Simulations. Macromol. Theory Simul. 2012, 21, 322−327. (19) Wisanrakkit, G.; Gillham, J. K. The Glass Transition Temperature as an Index of Chemical Conversion for a HighAmine/Epoxy System: Chemical and Diffusion-Controlled Reaction Kinetics. J. Appl. Polym. Sci. 1990, 41, 2885−2929. (20) Simpson, J. O.; Bidstrup, S. A. Rheological and Dielectric Changes During Isothermal Epoxy-Amine Cure. J. Polym. Sci., Part B: Polym. Phys. 1995, 33, 55−62. (21) Min, B. G.; Stachurski, Z. H.; Hodgkin, J. H. Cure Kinetics of Elementary Reactions of a Diglycidyl Ether of Bisphenol A/ Diaminodiphenylsulfone Epoxy Resin: 2. Conversion Versus Time. Polymer 1993, 34, 4488−4495. (22) Enns, J. B.; Gillham, J. K. Time−Temperature−Transformation (TTT) Cure Diagram: Modeling the Cure Behavior of Thermosets. J. Appl. Polym. Sci. 1983, 28, 2567−2591. (23) Min, B. G.; Stachurski, Z. H.; Hodgkin, J. H.; Heath, G. R. Quantitative Analysis of the Cure Reaction of DGEBA/DDS Epoxy Resins without and with Thermoplastic Polysulfone Modifier Using Near Infra-red Spectroscopy. Polymer 1993, 34, 3620−3627. (24) Jordan, C.; Galy, J.; Pascault, J.-P. Measurement of the Extent of Reaction of an Epoxy−Cycloaliphatic Amine system and Influence of the Extent of Reaction on its Dynamic and Static Mechanical Properties. J. Appl. Polym. Sci. 1992, 46, 859−871. (25) Marks, M. J.; Snelgrove, R. V. Effect of Conversion on the Structure−Property Relationships of Amine-Cured Epoxy Thermosets. ACS Appl. Mater. Interfaces 2009, 1, 921−926.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b01303. Various details about the steps to prepare atomistic models of cross-linked epoxy (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*(K.S.K.) E-mail: [email protected]. *(F.R.P.Jr.) E-mail: [email protected]. ORCID

Ketan S. Khare: 0000-0002-5487-5553 Frederick R. Phelan Jr.: 0000-0001-8004-5281 Author Contributions

K.S.K. and F.R.P.Jr. planned the study reported in this paper. K.S.K. performed the simulations under the supervision of F.R.P.Jr. Both authors were involved in writing the final manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS K.S.K. was supported at Georgetown University (GU) by the U.S. National Institute of Standards and Technology (NIST) Grant 70NANB14H214. K.S.K. is a Postdoctoral Fellow at GU and a Guest Researcher at NIST. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by U.S. National Science Foundation Grant ACI-1053575.66 Certain commercial materials are identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by NIST, nor does it imply that the materials identified are necessarily the best available for the purpose.

■

REFERENCES

(1) Materials Genome Initiative - Strategic Plan; Office of Science and Technology Policy: Washington, DC, 2014. (2) Pham, H. Q.; Marks, M. J. Epoxy Resins. In Ullmann’s Encyclopedia of Industrial Chemistry; Wiley-VCH Verlag GmbH & Co. KGaA: 2000. (3) Winey, K. I.; Vaia, R. A. Polymer Nanocomposites. MRS Bull. 2007, 32, 314−322. (4) Vaia, R. A.; Giannelis, E. P. Polymer Nanocomposites: Status and Opportunities. MRS Bull. 2001, 26, 394−401. (5) Baur, J.; Silverman, E. Challenges and Opportunities in Multifunctional Nanocomposite Structures for Aerospace Applications. MRS Bull. 2007, 32, 328−334. (6) Hogg, P. J. Composites in Armor. Science 2006, 314, 1100. J


Article

Macromolecules (26) Lin, P.-H.; Khare, R. Molecular Simulation of Cross-Linked Epoxy and Epoxy−POSS Nanocomposite. Macromolecules 2009, 42, 4319−4327. (27) Soni, N. J.; Lin, P.-H.; Khare, R. Effect of Cross-linker Length on the Thermal and Volumetric Properties of Cross-linked Epoxy Networks: A Molecular Simulation Study. Polymer 2012, 53, 1015− 1019. (28) Sirk, T. W.; Khare, K. S.; Karim, M.; Lenhart, J. L.; Andzelm, J. W.; McKenna, G. B.; Khare, R. High Strain Rate Mechanical Properties of a Cross-linked Epoxy Across the Glass Transition. Polymer 2013, 54, 7048−7057. (29) Soldera, A.; Metatla, N. Glass Transition of Polymers: Atomistic Simulation versus Experiments. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2006, 74, 061803. (30) Ferry, J. Viscoelastic Properties of Polymers, 3rd ed.; Wiley: New York, 1980; 641 pp. (31) Williams, M. L.; Landel, R. F.; Ferry, J. D. The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids. J. Am. Chem. Soc. 1955, 77, 3701−3707. (32) Bero, C. A.; Plazek, D. J. Volume-Dependent Rate Processess in an Epoxy Resin. J. Polym. Sci., Part B: Polym. Phys. 1991, 29, 39−47. (33) Lemay, J. D.; Swetlin, B. J.; Kelley, F. N. Structure and Fracture of Highly Cross-linked Networks. In Characterization of Highly Crosslinked Polymers; American Chemical Society: 1984; Vol. 243, pp 165− 183. (34) Choy, I.-C.; Plazek, D. J. The Physical Properties of BisphenolA-Based Epoxy Resins During and After Curing. J. Polym. Sci., Part B: Polym. Phys. 1986, 24, 1303−1320. (35) Plazek, D. J.; Choy, I. C. The Physical Properties of BisphenolA-Based Epoxy Resins During and After Curing. II. Creep Behavior Above and Below the Glass Transition Temperature. J. Polym. Sci., Part B: Polym. Phys. 1989, 27, 307−324. (36) Khare, R.; Paulaitis, M. E.; Lustig, S. R. Generation of Glass Structures for Molecular SImulations of Polymers Containing Large Monomer Units: Application to Polystyrene. Macromolecules 1993, 26, 7203−7209. (37) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671−680. (38) While it may be noted that the simulation temperatures are well above the degradation temperature for these materials, the use of these high temperatures allows us to accelerate processes and capture them at simulation accessible time scales. Such a tactic is possible because the mechanisms underlying degradation are not included in the classical models used in this work. (39) Registered Trademark of Hexion Inc. Corporation, 180 East Broad St., Columbus OH 43215. (40) As we note in section 1 of the Supporting Information, we have used the so-called active form of the Epon 1001F molecule (Figure S1) instead of Epon 1001F shown in Figure 1a in order to minimize the topological modifications needed to transform the model from the unreacted mixture to the cross-linked network. Such a substition has no effect on the chemistry of the cross-linked network. (41) Lin, P.-H.; Khare, R. Local Chain Dynamics and Dynamic Heterogeneity in Cross-Linked Epoxy in the Vicinity of Glass Transition. Macromolecules 2010, 43, 6505−6510. (42) Lin, P.-H.; Khare, R. Glass Transition and Structural Properties of Glycidyloxypropyl-Heptaphenyl Polyhedral Oligomeric Silsesquioxane-Epoxy Nanocomposites. J. Therm. Anal. Calorim. 2010, 102, 461− 467. (43) Khabaz, F.; Khare, K. S.; Khare, R. Temperature Dependence of Creep Compliance of Highly Cross-linked Epoxy. AIP Conf. Proc. 2014, 1599, 262−265. (44) Sirk, T. W.; Karim, M.; Khare, K. S.; Lenhart, J. L.; Andzelm, J. W.; Khare, R. Bi-modal Polymer Networks: Composition-Dependent Trends in Thermal, Volumetric and Structural Properties from Molecular Dynamics Simulation. Polymer 2015, 58, 199−208. (45) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174.

(46) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (47) General Amber Force Field (Version 1.8, March 2015) included in AmberTools 16. (48) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132−146. (49) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (50) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: A New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (51) Case, D. A.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Luchko, T.; Luo, R.; Madej, B.; Mermelstein, D.; Merz, K. M.; Monard, G.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; Botello-Smith, W. M.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Xiao, L.; Kollman, P. A. AMBER 2016, University of California, San Francisco, 2016. (52) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (53) Sun, H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338−7364. (54) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Institute of Physics Publishing: Philadelphia, 1988; pp 267− 301. (55) Andersen, H. C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34. (56) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (57) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (58) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (59) Shinoda, W.; Shiga, M.; Mikami, M. Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation Under Constant Stress. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 134103. (60) Patrone, P. N.; Dienstfrey, A.; Browning, A. R.; Tucker, S.; Christensen, S. Uncertainty Quantification in Molecular Dynamics Studies of the Glass Transition Temperature. Polymer 2016, 87, 246− 259. (61) Rohatgi, A. WebPlotDigitizer, 3.10; Austin, TX, 2016. (62) Jójárt, B.; Martinek, T. A. Performance of the General AMBER Force Field in Modeling Aqueous POPC Membrane Bilayers. J. Comput. Chem. 2007, 28, 2051−2058. (63) Sprenger, K. G.; Jaeger, V. W.; Pfaendtner, J. The General AMBER Force Field (GAFF) Can Accurately Predict Thermodynamic and Transport Properties of Many Ionic Liquids. J. Phys. Chem. B 2015, 119, 5882−5895. (64) Mobley, D. L.; Dumont, É.; Chodera, J. D.; Dill, K. A. Comparison of Charge Models for Fixed-Charge Force Fields: SmallMolecule Hydration Free Energies in Explicit Solvent. J. Phys. Chem. B 2007, 111, 2242−2254. (65) Cerrada, M. L.; McKenna, G. B. Physical Aging of Amorphous PEN: Isothermal, Isochronal and Isostructural Results. Macromolecules 2000, 33, 3065−3076. (66) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; Roskies, R.; Scott, J. R.; Wilkins-Diehr, N. XSEDE: Accelerating Scientific Data. Comput. Sci. Eng. 2014, 16, 62−74. K


Quantitative Comparison of Atomistic Simulations with Experiment for

Recommend Documents