Development of Force Field Parameters for Molecular Simulation of

ARTICLE pubs.acs.org/JCTC

Development of Force Field Parameters for Molecular Simulation of Polylactide James H. McAliley and David A. Bruce* Department of Chemical and Biomolecular Engineering, Clemson University, Clemson, South Carolina 29634-0909, United States

bS Supporting Information ABSTRACT: Polylactide is a biodegradable polymer that is widely used for biomedical applications, and it is a replacement for some petroleum-based polymers in applications that range from packaging to carpeting. Efforts to characterize and further enhance polylactide-based systems using molecular simulations have to this point been hindered by the lack of accurate atomistic models for the polymer. Thus, we present force field parameters specifically suited for molecular modeling of PLA. The model, which we refer to as PLAFF3, is based on a combination of the OPLS and CHARMM force fields, with modifications to bonded and nonbonded parameters. Dihedral angle parameters were adjusted to reproduce DFT data using newly developed CMAP dihedral cross terms, and the model was further adjusted to reproduce experimentally resolved crystal structure conformations, melt density, volume expansivity, and the glass transition temperature of PLA. We recommend the use of PLAFF3 in modeling PLA in its crystalline or amorphous states and have provided the necessary input files required for the publicly available molecular dynamics code GROMACS.

’ INTRODUCTION Polylactide, also called polylactic acid (PLA), is an important polymer for biomedical applications, because it is compatible with living cells and is biodegradable.1,2 Further, PLA is of interest as a commodity polymer, and is used especially in single-use packaging applications.3 PLA is an α-polyester, and the primary structure of its repeat unit is shown in Figure 1. Classical molecular force fields, such as CHARMM4 and OPLS,5 have been widely used in recent decades for simulating organic molecules, by and large with good success. However, neither force field has been parametrized specifically for the dihedral angles present in α-polyesters such as PLA. In particular, dihedral interaction parameters for the OSCCαOS, CCαOSC, and CαOSCCα motifs, all of which are unique to α-polyesters, are not found in the parameter databases for these force fields. In lieu of these specific four-atom interaction parameters, one would typically use the so-called wildcard parameters included in the force field (these are general parameters represented by XCCαY, XCαOSY, and XOSCY, where X and Y may be any atom type). Though wildcard parameters provide a reasonable guess for the dihedral interactions in cases where more accurate parameters are unavailable, it has been shown that use of the wildcard parameters for α-polyesters results in poor accuracy when modeling PLA because the wildcard rotational energy barriers are centered at the wrong dihedral angles and do not describe the barrier heights predicted via quantum (DFT and MP2 level) models well.6,7 In this work, we develop a classical force field model specifically suited for polylactides, based on the OPLS and CHARMM forms. The present force field follows the work of O’Brien,7 in which the PLAFF model was developed and validated extensively for crystalline PLA, and of McAliley,6 in which the model was further developed for accuracy in modeling amorphous PLA (the resulting force field was referred to as PLAFF2). Our present r 2011 American Chemical Society

model is PLAFF3, and it differs from previous versions in its use of the CMAP cross-term dihedral potential originally developed for CHARMM.8 This potential function provides more flexibility in fitting barriers to bond rotation and has allowed us to fit the glass transition temperature of PLA with far more accuracy than could be obtained from linear combinations of individual dihedral potentials. We demonstrate that PLAFF3 is better suited than its predecessors for modeling the amorphous dihedral angle distributions in PLA, and that it retains the accuracy of PLAFF in simulating the crystalline phase. We believe that these parameters will be of value to the biological science community in studying PLA. Further, with the growing interest in using renewable polymers for commodity packaging applications, this model will likely be of use to the materials science community in exploring new PLA-based materials. Though other work on PLA force field development has appeared in the literature,912 to our knowledge, the PLAFF3 parameters represent the first noncommercial molecular model validated against electronic structure calculations and experimental data for PLA in melt, glassy, and crystalline phases. As such, we hope this work will allow a larger number of researchers to study the material through simulation than was previously feasible.

’ METHODS The fitting procedure used in this work is shown schematically in Figure 2. The procedure begins with assembling target data and providing an initial guess for the force field parameters. As a first step, the torsional parameters are adjusted to match DFT data obtained in previous work.13 Next, the model is tested against experimental crystal structure data for PLA. Dihedral parameters are then adjusted accordingly, until reasonable Received: April 14, 2011 Published: September 19, 2011 3756

dx.doi.org/10.1021/ct200251x | J. Chem. Theory Comput. 2011, 7, 3756–3767

Journal of Chemical Theory and Computation agreement is obtained with the experimental data for crystalline PLA. Following this step, the model is used to simulate the polymer in its melt state. The volume expansivity, β, is estimated from these simulations and compared with experimental dilatometric measurements. In some cases, an adjustment of the relative energies between energy minima can affect a change in β, by establishing a different temperature dependence of the polymer’s rotational isomeric states. However, if a more drastic change is required, the nonbonded parameters are adjusted for those atom types that are unique to α-polyesters, until the density and volume expansivity are near experimental values. After such adjustments, the entire fitting procedure must be repeated to ensure the agreement with DFT and that crystal structure data is maintained. Finally, the model is used in quench simulations, where the polymer is rapidly cooled from the melt

Figure 1. Chemical structure of PLA. Superscript labels on atoms are used for reference in the text. The three main chain bond rotations are labeled according to the convention for polypeptides.

ARTICLE

state into the glassy state. Using the WilliamsLandelFerry (WLF) equation, the resulting glass transition temperature, Tg, may be compared to experimental measurements. If necessary, the energy barriers are then adjusted for rotation about each main chain dihedral angle until agreement is reached with experimental Tg values. Each of these steps will be discussed in greater detail in the following sections. Initial Force Field Parameters. We considered two force fields as a starting point for the PLAFF models: The Optimized Potentials for Liquid Simulations5 (OPLS) and the force field from the Chemistry at Harvard Molecular Mechanics (CHARMM) package.4 The OPLS parameters were taken from the OPLS-AA parameter files as distributed with GROMACS14 version 3.3.3. For the CHARMM force field, parameters were taken from the CHARMM27 proteinlipid parameter files distributed with CHARMM version c32b2. Atom types were assigned on the basis of chemical functionality (see Supporting Information). Partial atomic charges were unaltered in each force field, with the exception of main-chain atoms and the carbonyl oxygen, which were adjusted slightly to achieve charge neutrality in the lactyl residue and to improve agreement with DFT results. The needed CHARMM27 parameters for PLA were ported into GROMACS, and all further molecular mechanics calculations reported for the CHARMM force field were performed in GROMACS

Figure 2. Flow diagram showing the procedure for fitting PLA force field parameters. The dihedral angles for ϕ and ψ that define the gt and gc energy minima are shown later in Figure 4. 3757

dx.doi.org/10.1021/ct200251x |J. Chem. Theory Comput. 2011, 7, 3756–3767

Journal of Chemical Theory and Computation

ARTICLE

Figure 3. DFT potential energy profile for rotation of the OSC bond through all values of dihedral angle ω. Calculations performed on a PLA trimer in vacuo.13

version 3.3.3. It is important to note that the parameters, especially the nonbonded parameters, from the OPLS and CHARMM force fields are not fully compatible. This mixed set of parameters was only used as a starting point for further parameter optimization. In the fitting procedure, bond stretching and angle bending parameters were optimized to more accurately reproduce energy potentials predicted using DFT (see the Supporting Information). In general, DFT methods are known to give accurate geometries, whereas they are less accurate at predicting vibrational frequencies. For this reason, geometric parameters for bonds and angles (the b0 parameter for bonds and the θ0 parameter for angles) were fit to DFT data, but the bond and angle force constants (kb and kθ, respectively) were unaltered from their original OPLS values. In this way, we deviated as little as possible from the OPLS model. Three dihedral interactions were also adjusted to achieve better agreement with the bond rotational potential energy surfaces calculated from DFT. These correspond to the backbone dihedrals labeled as ϕ, ψ, and ω in Figure 1 and are defined by the IUPAC convention using the main chain atoms (OS, C, and Cα). In PLAFF3, the potential energies for these three dihedral interactions were represented by tabulated functions. For rotation about ω, a one-dimensional tabulated function was used (see the GROMACS User Manual14 for more information), whereas a two-dimensional tabulated function (also called a correction map or CMAP15) was used for each pair of neighboring ϕ, ψ dihedrals. The CMAP potential was recently implemented in GROMACS,8 and all calculations involving such terms were performed with GROMACS version 4.5.1. Target Data. We used several criteria to select target data for parameter fitting. One criterion was that the model should be consistent with results from higher-level molecular simulation methods, namely, our DFT results from previous work.13 In addition, we aimed to be consistent with experimental results. Because PLA is often used in its semicrystalline form, we desired a model that could reproduce the properties of both the crystalline and amorphous states of the material. Thus, conformational data for the crystalline form of PLA were used, as well as the glass transition temperature and volumetric data for the amorphous polymer. DFT Data. The target potential energy values were taken from in vacuo DFT calculations (B3LYP/6-31G**) for a methyl

Figure 4. (a) DFT potential energy surface, with dihedral angles ϕ and ψ as independent variables. Calculations performed on a PLA trimer in vacuo.13 The two lowest energy minima, gauche/trans (gt) and gauche/cis (gc) are indicated. (b) Molecular geometry of gc conformation. (c) gt conformation.

terminated PLA trimer (CH3(OC(O)CH(CH3))3OCH3), as reported previously.13 These include an estimate of the potential energy barriers encountered with bond stretching and angle bending (see the Supporting Information), during rotation of the OSC bond described by ω (reproduced in Figure 3), and likewise for simultaneous rotation of bonds described by ϕ and ψ (Figure 4) for the central repeat unit in a PLA trimer. Crystal Structure Data. Several studies on the crystal structure of PLA have appeared in the literature.1620 We have chosen to use the structural coordinates from Sasaki and Asakura20 as our target data. The authors’ use of the linked atom full-matrix leastsquares (LAFLS) method21 and the Rietveld whole-fitting method22 allowed for the positions of individual atoms in the unit cell to be determined with a high degree of accuracy. The authors derived the α-form of the crystal structure from WAXD data, resulting in a frustrated 103 helix. The orthorhombic (α = β = γ = 90°) unit cell from that study has P212121 symmetry and lattice constants of a = 10.66(1) Å, b = 6.16(1) Å, and c = 28.88(2) Å.20 Volumetric Data for Amorphous/Melt PLA. To represent the volumetric properties of amorphous PLA, we selected the experimental data from Sato et al.,23 where the specific volumes of polylactide samples were measured at various temperatures and pressures using metal bellows dilatometry. While the experimental data cover a wide range of temperatures, the only specific volume data used in this study are plotted in Figure 5, and these data correspond to specific volumes measured by heating PLA samples at 1 bar above the melting temperature, Tm. It can be seen in Figure 5 that an abrupt change in volume occurs upon heating above Tm, which is attributed to the change in volume that occurs when the crystallites in the sample become amorphous. In practice, we use molecular simulations to study the melt phase of polymers at temperatures higher than those shown in Figure 5, utilizing the well-known timetemperature superposition principle for polymers.24 Thus, we look at the volume expansivity (see eq 1 below) to facilitate a comparison. From the data 3758



ARTICLE

Table 1. Some Reported Values of the Glass Transition Temperature of PLA

Figure 5. Target volumetemperature data at 1 bar, taken from Sato et al.23 Arrows indicate the authors’ estimate of the glass transition temperature, Tg, and the melting temperature, Tm, taken from separate calorimetry data.

in Figure 5, an expansivity of β = 7.8 104 K1 is calculated for the melt state. Glass Transition Data. It is generally accepted that barriers to bond rotation play a major role in determining the glass transition temperature of polymers. Experimentally, the value of the PLA glass transition is dependent on the method used to measure it, and it has also been reported to vary widely with moisture content.25 Common reported values of Tg for PLA are in the range of 327 to 345 K (see Table 1), which were obtained using differential scanning calorimetry (DSC) and dielectric relaxation spectroscopy (DRS). The value reported by Auras was measured after extensive drying of the PLA samples.25 Since water is known to have a plasticizing effect on the material, it follows that this estimate is at the high end of the reported range of Tg values. Fitting Procedure Using DFT Target Data. In fitting our model to DFT data, we performed energy minimizations using the force field model, each constrained by the same independent variable(s) as the DFT energy minimizations (i.e., ω or ϕ and ψ) before comparing energies with the DFT results. This energy minimization step with the force field model adds a high level of nonlinearity to the fitting procedure. Each time the dihedral parameters are adjusted, the minimum energy conformation at each independent variable also changes. Thus, obtaining the optimal torsional potentials according to this prescription requires an iterative scheme. We began each iteration of the fitting procedure with the dihedral having the largest potential energy barriers. Thus, the rotational energy barrier for the OSC bond (dihedral angle ω) was fit to the DFT data shown in Figure 3 first, followed by a simultaneous fit of the ϕ and ψ potentials to the data shown in Figure 4. A weighted least-squares approach was used in developing PLAFF2, as described elsewhere.6 However, in PLAFF3, simple adjustments to the tabulated potentials were used to match the DFT data. Weighting of data points was not necessary due to the greater flexibility of the tabulated potentials versus the traditional cosine expansions used in PLAFF2. After each parameter optimization step, the force field minimum energy conformations were reevaluated using the most current dihedral parameters. This process was repeated until reasonable convergence was achieved with respect to the dihedral parameters and the minimized energies.

lead author

method

rate

Tg (K)

26

Dorgan

DSC

10 °C/min

331.6

Sato23 Auras25

DSC DSC

10 °C/min

337 344.6

Joziasse27

DSC

10 °C/min

336

Kanchanasopa28

DRS

τ = 100 s

327

Bounded Adjustment of Dihedral Potentials. We employ a bounded adjustment procedure for altering the dihedral potentials, to avoid drastically changing the force field parameters and thus the minimum energy conformations. For all iterations, the tabulated potentials were adjusted to be as close as possible to the DFT target data without exceeding a specified change in energy. Otherwise, if unbounded adjustments were allowed, we often observed divergent behavior due to the nonlinear aspects of the iterative scheme. It was found that suitable stability was achieved by limiting the change in energy at each tabulated point to 10 kJ/mol for fitting the ω dihedral, and 5 kJ/mol for fitting ϕ and ψ dihedrals. Because these limits gave satisfactory performance, no attempt was made at further tuning the fitting procedure with respect to them. Refinement Using Crystal Structure Data. Returning to our discussion of Figure 2, we proceed, after sufficient convergence is obtained in fitting to DFT bond rotation data, by examining the crystal structure of PLA with the resulting force field parameters. In these simulations, a super cell based on the crystalline unit cell is built according to the WAXD-resolved structure of Sasaki and Asakura.20 This super cell contains 32 PLA chains, each containing 50 monomers, with the boundary conditions for the polymer being such that monomer number 1 in a chain was bonded to monomer number 50 from the neighboring periodic cell. The selected system size is sufficiently large that no finite size effects are observed with the simulated systems. The system is simulated for 3.0 ns in the NPT ensemble, whereby the lattice or box dimensions are allowed to adjust to their equilibrium values; however, all lattice angles (α, β, and γ) were constrained to 90°, so as to maintain an orthorhombic unit cell that matched experimental observations. Anisotropic pressure coupling was applied with the Berendsen algorithm, such that each box length was adjusted independently.29 The NoseHoover thermostat was used to control temperature at 300 K.30,31 A cutoff of 1.0 nm was used for van der Waals interactions, while the electrostatics were treated with the Particle-Mesh Ewald (PME) method.32 While one of our goals was to have a minimized PLA system that accurately reproduced the experimentally observed dihedral angles from diffraction studies, obtaining these in the crystal structure is difficult without first having accurate bond lengths and valence angles due to packing considerations in the unit cell. Therefore, before earnestly examining the predicted dihedral values, we optimize the bond and angle force field parameters using a series of in vacuo DFT calculations that examine the variation in PLA system energy as a function of perturbations to each bond and angle from its minimum energy value (see Supporting Information). This is done in much the same way as fitting the dihedral parameters to DFT bond rotation data, using a self-consistent iteration scheme. Once the prerequisite of accurate valence geometries is achieved, the dihedral angles (ϕ, ψ, and ω) and lattice vectors 3759


Journal of Chemical Theory and Computation are examined over the final 1.0 ns of dynamics of a crystal structure simulation, and their values are compared to those reported in the experimental literature. Should the simulations be inconsistent with the experimental data, the position of the gt energy minimum (see Figure 4) is adjusted with respect to ϕ and ψ, as is the position of the trans energy minimum for ω. This process is repeated until the experimental dihedral values are accurately reproduced. Refinement Using Melt Phase Target Data. Once adequate agreement with the crystalline unit cell is obtained, simulations are carried out on amorphous PLA using isothermalisobaric replica exchange molecular dynamics (NPT-REMD) as implemented in GROMACS. Four separate NPT-REMD runs were performed with unique input configurations, which were initially generated using the Amorphous Cell module in Accelrys’ Materials Studio version 4.4.6 In our implementation, each replica is comprised of three chains, each containing 500 repeat units (refer to Figure 1), and two lactide molecules. The chain length was chosen to be greater than the experimental entanglement length, which is approximately 125 repeat units. Additionally, limited simulations with other chain length PLA systems showed that the calculated Tg value decreased with chain length as expected, but no detailed investigation of this effect was attempted. Lactide molecules were also included in the simulated systems because there is always a small percentage of residual lactide monomer in real polylactide samples, and these monomers have a plasticizing effect on the material. With two lactide molecules per simulation cell, our simulated PLA system contains 0.26% residual lactide on a weight basis; the specific amount of lactide present in an industrially produced PLA resin is usually less than 1%,33 and 0.2 to 0.3 wt %34 is common. In each replica exchange simulation, the average volume is calculated for the melt state as a function of temperature. From this, the volume expansivity can be estimated graphically using the relation35 1 ∂v β¼ ð1Þ v ∂T P where v is the specific volume of the system. Since each replica has the same pressure, a plot can be constructed of ln v versus T, and β can be estimated from its slope. This is compared to experimental measurements of the expansivity of PLA. If satisfactory agreement is not obtained, this indicates that the nonbonded parameters may need further adjustment. When it was necessary to alter the nonbonded parameters, the atom types for the ester oxygen (OS) and α-carbon (Cα) were chosen for adjustment. These atom types were selected because they are the most likely to deviate from the behavior of normal esters, and no such atom types exist in OPLS or CHARMM for α-polyesters. Note that, when adjustment of these atoms’ nonbonded parameters is necessary, the dihedral angles must again be readjusted to preserve agreement with the crystal structure. Refinement Using Glass Transition Target Data. The glass transition temperature, Tg, is commonly interpreted for polymers as the temperature below which bond rotations are kinetically trapped. That is, it is the temperature below which torsional energy barriers are crossed at rates much longer than the time scale on which the polymer is observed. As such, the value of Tg for PLA is influenced by the height of the energy barrier between the various rotational isomeric states.

ARTICLE

Many studies have appeared in the literature examining the glass transition temperature via molecular dynamics,3638 though relatively few papers address the temporal dependence of the observed glass transition temperature.39,40 It is well-known, experimentally, that the glass transition will be observed at higher temperatures when a polymer specimen is cooled at a faster rate.41 This behavior is described very well, over the range of experimental time scales, by the WLF equation:42 B ðT T0 Þ f0 ð2Þ ln AT ¼ f0 þ ðT T0 Þ αf where AT is the reduced variables shift factor, B is a constant, f0 is the fractional free volume of the polymer at the reference temperature T0, and αf is the coefficient of expansion of the free volume. Although the quenching rates accessible to molecular dynamics simulations can differ from experimental cooling rates by 14 orders of magnitude or more, the validity of the WLF equation over such wide a temporal range has been established recently through molecular simulation.40 The glass transition temperature of PLA was estimated from our force field model by quenching the amorphous conformations from the NPT-REMD simulations, using a replica at 604.5 K as the starting structure. Simulation conditions were identical to those used in each of the NPT-REMD replicas, except that the set point of the Nose-Hoover thermostat was varied linearly with simulation time over the entire run. Each run lasted until a temperature of 300 K was reached. From each of the four NPTREMD simulations, six separate quench runs were performed, with quench rates of 15 K/ns, 30 K/ns, 60 K/ns, 150 K/ns, 300 K/ns, and 600 K/ns. The glass transition temperature was estimated for each run by fitting a straight line to a plot of ln v versus T, using all data points below 400 K. A second straight line was drawn through the melt data taken from the NPT-REMD runs, for temperatures above 500 K. The intersection of the two lines was taken as Tg. Such estimates were then averaged for each quench rate, and then a least-squares fit was performed using the WLF model (eq 2). The reference quench rate was taken to be 10 K/min (normal lab conditions for measuring Tg). We found that the so-called universal WLF constants (B/f0 = 40.16 and f0/αf = 51.6 K) provided a good fit with the PLA force-field models, reducing the WLF equation to a single adjustable parameter, T0. When used in this manner, T0 corresponds to the Tg observed at the reference quench rate.

’ RESULTS AND DISCUSSION In implementing the iterative procedure described in Figure 2, many intermediate sets of force field parameters were examined in this study. Here, we compare and discuss the important results from seven different models, ranging from the unaltered CHARMM and OPLS force fields to our intermediate parameter sets derived from those force fields and from the first generation PLAFF force field of O’Brien to the present version of our PLA force field, PLAFF3. Abbreviations for these different models are summarized in Table 2 for ease of reference. These include the models obtained directly after a least-squares fit to the DFT potential energies, referred to as the OPLS0 and CHARMM0 models. These two models demonstrate that fitting to the DFT energies alone is not sufficient to reproduce experimental data, 3760



ARTICLE

Table 2. Description of the Various Classical Models Discussed in the Texta model

description

OPLS

the OPLS force field as developed by Jorgensen and co-workers5,43 (all-atom version, also known as OPLS-AA)

CHARMM

the CHARMM force field as developed by Brooks and co-workers4

0

OPLS

OPLS, with backbone torsional potentials refit to DFT data

CHARMM0

CHARMM, with backbone torsional potentials refit to DFT data

OPLS00

OPLS, with CHARMM nonbonded parameters substituted for O1 and C3, selected bond stretching and angle bending terms

PLAFF

The PLA force field developed by O’Brien7 (all-atom version, also known as PLAFF-AA)

PLAFF2

the second version of PLAFF, with improved accuracy for simulating amorphous PLA;6 the model is OPLS00 , with backbone torsional parameters further adjusted to reproduce crystal structure data and

PLAFF3

the PLA force field developed in this work; the model is OPLS00 , with backbone

refit to DFT, and backbone torsional potentials refit to DFT data

to improve agreement with the experimental glass transition temperature of PLA torsional parameters further adjusted using the CMAP potential to reproduce crystal structure data and to improve agreement with the experimental glass transition temperature of PLA a

PLAFF3 is the recommended potential for molecular simulation of PLA.

and further adjustment was required as described in Figure 2. The OPLS00 model shows results in which the parameters have been further adjusted to match the bond lengths and valence angles reported in crystal structure studies of PLA, and in which the nonbonded parameters were adjusted by trial-and-error to better match the melt density and volumetric expansivity of PLA. Finally, the PLAFF2 and PLAFF3 models improve upon OPLS00 by adjusting the dihedral parameters to match crystal structure data and glass transition temperature data. The principal difference between PLAFF2 and PLAFF3 is that PLAFF3 uses CMAP dihedral cross terms, whereas PLAFF2 uses linear combinations of single dihedral potentials. Comparison of the Classical Models to DFT Data. The energy landscapes for bond rotation about ϕ and ψ are shown in Figure 6, calculated using the various models described in Table 2. The figure also shows the DFT target data for comparison. While we do not expect the optimum force field to be in complete agreement with DFT data, we desire the overall shape and location of relative minima/maxima to coincide with the DFT results in order to realistically model the amorphous configuration distribution. Thus, as shown in Figure 6b, the first generation PLAFF force field raises some concern, due to the presence of a low-energy local minimum in the vicinity of (ϕ, ψ) = (30°, 150°) that does not appear in the DFT potential energy surface. Additionally, the gc minimum is predicted by PLAFF to be a much less probable configuration than predicted by DFT. The presence of the extra minimum in Figure 6b is only of concern for applications involving amorphous phases of PLA, in which case the entire dihedral space may be accessed by the simulated polymer chains according to the energetics of the force field model. While the nonphysical local minimum near (ϕ, ψ) = (30°, 150°) is a striking feature of Figure 6b, it is also obvious from the figure that O’Brien was very successful in fitting the potential energy surface in the vicinity of the global minimum (in the gt position shown in Figure 4). This is evidenced by the remarkable performance of PLAFF in simulating crystalline PLA,7 and therefore, we feel that the original PLAFF is still very well-suited in modeling the crystalline phase of PLA. When examining the OPLS and CHARMM models in Figure 6c and d, we see that both models lack adequate representation of the global gt minimum predicted by DFT. This observation helps to explain

the superior performance of PLAFF in the crystalline phase as compared with OPLS and CHARMM and suggests that OPLS and CHARMM should not be used for crystalline or amorphous phase simulations without first correcting the backbone torsional potentials. Figure 6e and f show the results of performing a least-squares fitting procedure to alter the torsional potentials of OPLS and CHARMM, while leaving all other interaction parameters in the models unchanged. This figure demonstrates that there are limitations inherent in each model, preventing a perfect fit to the desired potential energy surface. For example, the CHARMM0 potential energy surface in Figure 6f still shows remnants of the local minima, situated in the negative ϕ region between the g c and g t energy minima of the CHARMM model in Figure 6d. The major shortcoming of the models shown in Figure 6bh is that corrections to the (ϕ, ψ) potential energy surface are limited to linear combinations of separate functions of ϕ and functions of ψ. Without the use of more sophisticated potential energy functions, e.g., the CMAP dihedraldihedral cross terms available in recent versions of the CHARMM program,8,15 accurately reproducing the entire two-dimensional potential energy surface of Figure 6a is highly dependent on the other interactions within the model, such as the bond stretching and angle bending parameters. Figure 6g and h give the potential energy surfaces after fitting the model to crystal structure data and glass transition data, respectively (see discussion in the following sections). The last plot (Figure 6i) is our currently recommended model, PLAFF3. Note that the agreement between PLAFF2 and the DFT data is diminished when compared to OPLS0 , since adjustments to fit one set of target data inevitably alters the performance of the model in reproducing all other target data. The resulting model is a compromise between competing target data. However, the addition of the CMAP dihedral term in PLAFF3 allows for a much more accurate fit of the potential energy surface, which can be altered in very specific local regions without affecting the shape of the surface elsewhere. For example, the gt global minimum of PLAFF3 was shifted by approximately +10° and 15° in ϕ and ψ, respectively, to improve agreement with the crystal structure, yet the remaining portions of the surface are unaffected by this shift. 3761



ARTICLE

Figure 6. Bond rotational energy profiles for the ϕ and ψ dihedrals (shown in Figure 1) of a PLA trimer, calculated from (a) B3LYP/6-31G**,13 (b) PLAFF,7 (c) OPLS,5 (d) CHARMM,4 (e) OPLS0 , (f) CHARMM0 , (g) OPLS00 , (h) PLAFF2,6 and (i) PLAFF3. Refer to Table 2 for a description of the models.

Table 3. Lattice Dimensions of PLA at 300 K from Published Studies and from Crystal Structure Simulationsa c (Å)

diff (%)

density (g/cm3)

diff (%)

5

28.88 29.01

1

1.261 1.472

16.7

6.10

1

28.80

0

1.285

1.8

0

6.45

5

27.80

4

1.247

1.2

10.46

1.9

6.05

1.8

31.14

7.8

1.214

3.8

10.72

0.6

5.97

3.1

31.47

9.0

1.188

5.8

10.51

1.4

5.97

3.1

31.36

8.6

1.216

3.6

8.78

17.6

6.03

2.1

34.67

20.0

1.303

3.3

OPLS00 PLAFF2

10.54 10.59

1.1 0.7

6.08 6.25

1.3 1.5

30.85 29.74

6.8 3.0

1.210 1.215

4.1 3.7

PLAFF3

10.70

0.4

6.14

-0.3

29.98

3.8

1.214

3.7

a (Å)

diff (%)

b (Å)

Sasaki20 Aleman16

10.66 9.66

9

6.16 5.80

Hoogsteen18

10.60

1

de Santis17

10.70

OPLS CHARMM OPLS0 CHARMM0

diff (%)

a

Refer to Table 2 for a description of the models. Differences are calculated with respect to the experimental study of Sasaki and Asakura.20 The recommended model, PLAFF3, is emphasized in bold.

Comparison of the Classical Models to Crystal Structure Data. Results from crystal structure simulations using each of

the models are shown in Table 3 and in Figure 7. In each case, the simulation results are compared to reference values from the 3762



ARTICLE

Figure 7. Dihedral angle distributions for crystalline PLLA at 300 K, simulated with (a) OPLS,5 (b) CHARMM,4 (c) OPLS0 , (d) CHARMM0 , (e) PLAFF2, (f) PLAFF3. Refer to Table 2 for a description of the models. Vertical dotted lines, values from the WAXD crystal structure analysis of Sasaki and Asakura;20 vertical dashed lines, averaged values from the PLAFF simulations performed by O’Brien.7

experimentally resolved crystal structure(s). We have given priority to ensure PLAFF3 matches the data in Table 3, as the crystalline density and lattice vectors can be measured directly with very few assumptions involved in the experimental analysis. Thus, the dihedral angle distributions in Figure 7 were considered a secondary target. Figure 7 shows the dihedral angle distributions during simulation of crystalline PLA. For each model, five different histograms are accumulated for each backbone dihedral, ω, ϕ, and ψ. These separate histograms are presented for each of the five unique residues in the frustrated helical structure predicted by Sasaki and Asakura.20 In each model, it is evident that these five residues take on different dihedral values, according to their orientation inside the unit cell. This supports the existence of a frustrated structure and demonstrates that a helix with perfect screw symmetry is not possible under the crystalline packing conditions of PLA. From Figure 7a and b, it is apparent that the OPLS and CHARMM models do not predict the same dihedral angle distribution as suggested by the WAXD results.20 A more surprising result was that refitting the torsional potentials to DFT data had very little effect on the dihedral angle distributions in the crystalline phase, as evident in Figure 7c and d. We found it essential to improve agreement with the

experimental unit cell lattice vectors before adjusting the dihedral parameters. The unit cell dimensions impose constraints on the set of dihedral angles that are probable, given that the crystal structure must be periodic with respect to those dimensions. Further, the set of bond lengths and angles played a vital role in achieving agreement with the crystal structure, as these impose the same sort of constraints on the dihedral angles when a periodic cell is used. Adjustments to bond stretching and angle bending parameters for this purpose are provided in the Supporting Information. Once the bonded interactions were adjusted and more closely matched those used in the WAXD analysis of Sasaki and Asakura,20 adjustment of the dihedral angles in the crystalline structure was relatively simple; in practice, we found that all of the backbone dihedral angle distributions could be shifted toward the WAXD values, by altering the potential with respect to the ϕ and ψ dihedral angles alone. A simple shift in the position of the global minimum was required, as depicted in Figure 8. The minimum was shifted by 10° in the ϕ dihedral angle and 15° in the ψ dihedral angle. A similar shift in ϕ was also required by O’Brien in developing PLAFF.7 In developing the PLAFF3 model, one of our stated goals was to obtain a force field that is suitable for modeling PLA in its 3763



ARTICLE

Table 4. Volume Expansivities Estimated for Melt Phase PLAa method/model OPLS

β 104 (K1) 7.74 ( 0.07

CHARMM PLAFF3

9.5 ( 0.1 7.7 ( 0.4

Sato et al.23

7.8 ( 0.4

a

Values are calculated from the simulation results shown in Figure 9, by a linear regression (on a log scale plot) of the data points above 550 K. An estimate using the experimental data of Sato et al.23 is included for comparison. Listed errors are 95% confidence intervals for each slope.

Figure 8. Adjustment of the torsional potential for the ϕ and ψ dihedral angles, which resulted in improvement of the dihedral angle distributions in crystalline simulations.

Figure 10. Example of the specific volumetemperature (vT) plot used to determine the glass transition temperature. Results are from the OPLS model, using two different quench rates.

Figure 9. Melt phase densities of PLA, plotted from four separate NPTREMD simulations for each of the CHARMM, OPLS, and PLAFF3 models. The melt phase experimental measurements of Sato et al.23 are included for comparison and extrapolated toward the higher simulation temperatures.

amorphous state. Simultaneously, we wished to retain the model’s accuracy in simulating crystalline PLA, which was a hallmark of O’Brien’s original PLAFF.7 We believe the results presented thus far demonstrate that PLAFF3 does indeed accurately predict the crystalline structure of PLA. In addition, Figure 6i shows the improvement in the topography of the PLAFF3 bond rotational energy landscape, when compared to PLAFF and PLAFF2, and demonstrates that the new model is more likely to have the correct dihedral angle distribution in the melt and amorphous state. In what follows, we show that PLAFF3 is also better suited for simulating PLA in its noncrystalline form,

up to high temperatures when compared to the other models discussed here. Comparison of the Classical Models to Melt Phase Dilatometric Data. When examining the models’ performance in high temperature simulations, we found that the OPLS model underpredicts the specific volume of PLA in the melt phase. This is shown in Figure 9, using results from the NPT-REMD simulations. CHARMM, on the other hand, tends to overestimate the specific volume. Results from the OPLS-based force fields generally reproduced the volume expansivity of PLA, as shown in Table 4, whereas the CHARMM-based models tended to have higher expansivities than indicated in the experimental results of Sato et al.23 It was found that substituting one or more of the nonbonded parameters (both Lennard-Jones parameters and partial charges) from CHARMM helped to increase the specific volume in the melt, without increasing the expansivity above the desired range. Following this observation, in the PLAFF3 force field, CHARMM nonbonded parameters are used for the OS and Cα atoms. While still slightly lower than the experimental measurements, the melt volumes predicted by PLAFF3 are noticeably closer to the experiment than either OPLS or CHARMM; this result supports our assertion that the model may be used equally well in simulating the melt and/or crystalline states of PLA. Comparison of the Classical Models to Glass Transition Data. The last material property we used in constructing the PLAFF3 set of parameters was the PLA glass transition temperature, Tg. Figure 10 gives an example of the specific volume intersection method used for determining Tg, at two different quench rates, using the OPLS force field. The results depicted in 3764



ARTICLE

Figure 11. WLF plot for extrapolating simulation glass transition data to realistic (laboratory scale) quench rates. Results are from the PLAFF3 model. Universal WLF constants are used, with a lab-scale quench rate of q0 = 10 K/min. Error bars are propagated from 95% confidence intervals on the slopes and intercepts of the melt and glassy vT plots (see Figure 10). Here, T0 = 327 K is the glass transition temperature extrapolated to the lab-scale quench rate.

Table 5. Glass Transition Temperatures Calculated from the Various Models Explored in This Worka Tg (K)

method/model simulation data

388 ( 14

OPLS CHARMM

367 ( 15

OPLS00 PLAFF44

403 ( 12 408b

PLAFF26

386 ( 11

PLAFF3

327 ( 12 Experimental Data

Dorgan26

331.6

Kanchanasopa28

327

a

Simulation results from previous studies using PLAFF44 and PLAFF26 are included for reference, as well as selected experimental results25,28 for PLA. b Extrapolated to infinite molecular weight limit, not corrected for quench rate dependence.

the figure are generally representative of all intersection plots constructed during this work; the faster quenching rates consistently gave intersection points that are higher up on the melt volumetric curve. Figure 11 demonstrates the extrapolative method used to estimate Tg for laboratory scale quench rates using the WLF relation. In fitting the WLF equation to the simulation data in Figure 11, the only adjustable parameter used was T0, which corresponds to the laboratory-scale glass transition temperature when the universal WLF constants are used (see the Methods section of this paper). A survey of the glass transition temperatures for some of the models discussed in this work is presented in Table 5. Not all models were tested for the glass transition temperature; following our procedure laid out in Figure 2, we required that our models perform accurately in both the crystalline and melt states before attempting to examine the glass transition temperature. Thus, the OPLS0 and CHARMM0 models were not examined

with glass transition simulations, as they did not meet the prerequisites in simulating the crystal structure. Similarly, PLAFF was not used because it is believed to give inadequate dihedral angle distributions. We made three exceptions, for demonstration purposes. We chose to estimate Tg using OPLS, CHARMM, and OPLS00 , because these results give some idea of how the glass transition temperature was affected by changes made early on in the fitting procedure. Most of the simulation-based estimates of Tg shown in Table 5 are higher than the experimentally observed glass transition temperature, with the PLAFF3 force field being the closest to the experimental value. It is also apparent in Table 5 that the modification of the torsional and other potentials from the OPLS to the OPLS00 model resulted in a worsening of the Tg estimate using OPLS00 . It is obvious that, in adjusting the nonbonded and valence interactions in OPLS to obtain the OPLS00 model, we affected the barrier height of bond rotation about the ψ dihedral angle. In the PLAFF2 model, we were able to remove this artifact, yet the limitations of using uncorrelated (non-CMAP) dihedral potentials for ϕ and ψ made further lowering of the barriers difficult without drastically affecting the overall potential energy surface with respect to ϕ and ψ, as discussed elsewhere.6 By introducing the CMAP potential in PLAFF3, we were able to remedy this problem. The considerable freedom entailed in the CMAP model allowed for a nearly exact fit of the DFT potential energy surface in Figure 4, and any adjustments made during crystal structure fitting could be made independently of the bond rotational barrier heights. Whereas such changes in PLAFF2 resulted in an unwanted increase in the barrier to rotation about the ψ dihedral angle, this barrier height could be preserved using the CMAP potential in PLAFF3. We found that the bond rotational barrier height in PLAFF3 was very close to the DFT-predicted value. PLAFF3 predicted a glass transition temperature within the range of the experimental results, without need for further adjustment of the barrier to rotation about the ψ dihedral angle. This is a major improvement over previous versions of the force field.

’ CONCLUSIONS In this paper, we presented our work related to the development of an optimized model for the atomistic simulation of polylactide (PLA). The model, PLAFF3, was shown to perform well in simulations of the amorphous and crystalline states of PLA. This model is an update to the previous versions by O’Brien7 and McAliley,6 and we have significantly improved the ability of the model to describe the proper dihedral angle distributions in the amorphous states of PLA. On the basis of the results of this work, we recommend the use of the PLAFF3 model under most circumstances. A major improvement in PLAFF3 over previous models is its ability to predict the glass transition temperature of PLA. This was possible due to the CMAP dihedral cross terms that were used in PLAFF3, in place of the linear combinations of individual dihedral terms used in prior versions of the force field. The inability to reproduce the experimental glass transition temperature was the largest shortcoming of PLAFF2, and PLAFF3 addresses this problem while retaining the accuracy of PLAFF2 in simulating the melt and crystalline states of PLA. We believe the wide range of properties captured by the PLAFF3 model make it well suited for studying a wide range of phenomena, such as crystallization, permeant diffusion, and shear and elongational 3765


Journal of Chemical Theory and Computation flow. Due to the accuracy of PLAFF3 in simulating the pure crystalline form, the model should work equally as well as PLAFF in simulating various surface interactions with the polymer, since the polymer is known to exhibit a high level of crystallization at surface boundaries. We feel that a judicious practitioner of molecular modeling should be able to apply PLAFF3 to successfully simulate any of these phenomena on a molecular level.

’ ASSOCIATED CONTENT

bS

Supporting Information. Parameters for the PLAFF3 force field as well as complete simulation input files for PLA. The input files, including atomic coordinate files, are specifically for the GROMACS molecular dynamics package and employ the PLAFF3 force field; however, the parameters in these files can be easily ported to other simulation packages. The atomic coordinate files are for fully equilibrated melt configurations of PLA. Additionally, the DFT derived energy profiles associated with bond stretching and angle bending of an in vacuo PLA trimer are provided. This information is available free of charge via the Internet at http://pubs.acs.org

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT The authors acknowledge the support of the Center for Advanced Engineering Fibers and Films and the ERC program of the National Science Foundation (NSF) under Award Number EEC-9731680. Partial support was also provided by the National Institutes of Health (NIH)/National Institute of Bioimaging and Bioengineering (NIBIB) under grant number R01 EB006163. The authors would also like to acknowledge the support of the staff from the Cyberinfrastructure Technology Integration group, and the use of the advanced computational resources deployed and maintained by Clemson Computing and Information Technology. ’ REFERENCES (1) Middleton, J. C.; Tipton, A. J. Synthetic biodegradable polymers as orthopedic devices. Biomaterials 2000, 21 (23), 2335–2346. (2) Anderson, J. M.; Shive, M. S. Biodegradation and biocompatibility of PLA and PLGA microspheres. Adv. Drug Delivery Rev. 1997, 28 (1), 5–24. (3) Auras, R; Harte, B; Selke, S. An overview of polylactides as packaging materials. Macromol. Biosci. 2004, 4 (9), 835–864. (4) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. Charmm - a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4 (2), 187–217. (5) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225–11236. (6) McAliley, J. H. Development of improved torsional potentials in classical force field descriptions of poly (lactic acid). Ph.D. Dissertation, Clemson University, Clemson, SC, 2009. http://etd.lib.clemson.edu/ documents/1252938067/ (accessed Dec 2010) (7) O’Brien, C. P. Quantum and molecular modeling of polylactide. Ph.D. Dissertation, Clemson University, Clemson, SC, 2005.

ARTICLE

(8) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comp. 2010, 6 (2), 459–466. (9) Blomqvist, J. Ris metropolis monte carlo studies of poly(l-lactic), poly(l,d-lactic) and polyglycolic acids. Polymer 2001, 42 (8), 3515–3521. (10) Blomqvist, J.; Ahjopalo, L.; Mannfors, B.; Pietila, L. O. Studies on aliphatic polyesters i: Ab initio, density functional and force field studies of esters with one carboxyl group. THEOCHEM 1999, 488, 247–262. (11) Blomqvist, J.; Mannfors, B.; Pietila, L. O. Studies on aliphatic polyesters. Part ii. Ab initio, density functional and force field studies of model molecules with two carboxyl groups. THEOCHEM 2000, 531, 359–374. (12) Blomqvist, J.; Mannfors, B.; Pietila, L. O. Amorphous cell studies of polyglycolic, poly(l-lactic), poly(l,d-lactic) and poly(glycolic/l-lactic) acids. Polymer 2002, 43 (17), 4571–4583. (13) McAliley, J. H.; O’Brien, C. P.; Bruce, D. A. Continuum electrostatics for electronic structure calculations in bulk amorphous polymers: Application to polylactide. J. Phys. Chem. A 2008, 112 (31), 7244–7249. (14) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible and Free. J. Comput. Chem. 2005, 26, 1701–1718. (15) Mackerell, A. D.; Feig, M.; Brooks, C. L. Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25 (11), 1400–1415. (16) Aleman, C.; Lotz, B.; Puiggali, J. Crystal structure of the alphaform of poly(l-lactide). Macromolecules 2001, 34 (14), 4795–4801. (17) de Santis, P.; Kovacs, J. Molecular conformation of poly(s-lactic acid). Biopolymers 1968, 6 (3), 299–306. (18) Hoogsteen, W.; Postema, A. R.; Pennings, A. J.; Tenbrinke, G.; Zugenmaier, P. Crystal-structure, conformation, and morphology of solution-spun poly(l-lactide) fibers. Macromolecules 1990, 23 (2), 634–642. (19) Kobayashi, J.; Asahi, T.; Ichiki, M.; Oikawa, A.; Suzuki, H.; Watanabe, T.; Fukada, E.; Shikinami, Y. Structural and optical-properties of poly lactic acids. J. Appl. Phys. 1995, 77 (7), 2957–2973. (20) Sasaki, S.; Asakura, T. Helix distortion and crystal structure of the alpha-form of poly(l-lactide). Macromolecules 2003, 36 (22), 8385–8390. (21) Arnott, S.; Wonacott, A. J. Atomic coordinates for an alphahelix refinement of crystal structure of alpha-poly-l-alanine. J. Mol. Biol. 1966, 21 (2), 371–383. (22) Rietveld, H. M. A profile refinement method for nuclear and magnetic structures. J. Appl. Crystallogr. 1969, 2 (2), 65–71. (23) Sato, Y.; Inohara, K.; Takishima, S.; Masuoka, H.; Imaizumi, M.; Yamamoto, H.; Takasugi, M. Pressure-volume-temperature behavior of polylactide, poly(butylene succinate), and poly(butylene succinateco-adipate). Polym. Eng. Sci. 2000, 40 (12), 2602–2609. (24) Painter, P. C.; Coleman, M. M. Essentials of polymer science and engineering; DEStech Publications: Lancaster, PA, 2009; p 464. (25) Auras, R.; Harte, B.; Selke, S. An overview of polylactides as packaging materials. Macromol. Biosci. 2004, 4 (9), 835–864. (26) Dorgan, J. R.; Lehermeier, H.; Mang, M. Thermal and rheological properties of commercial-grade poly(lactic acid)s. J. Polym. Environ. 2000, 8 (1), 1–9. (27) Joziasse, C. A. P.; Veenstra, H.; Grijpma, D. W.; Pennings, A. J. On the chain stiffness of poly(lactide)s. Macromol. Chem. Phys. 1996, 197 (7), 2219–2229. (28) Kanchanasopa, M.; Runt, J. Broadband dielectric investigation of amorphous and semicrystalline l-lactide/meso-lactide copolymers. Macromolecules 2004, 37 (3), 863–871. (29) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684–3690. (30) Nose, S. A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81 (1), 511–519. (31) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31 (3), 1695–1697. 3766



ARTICLE

(32) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - An N 3 Log(N) Method For Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089–10092. (33) Kolstad, J. J.; Witzke, D. R.; Hartmann, M. H.; Hall, E. S.; Nangeroni, J. Lactic Acid Residue Containing Polymer Composition and Product Having Improved Stability, and Method for Preparation and Use Thereof. U.S. Patent 6,353,086, March 5, 2002. (34) Bopp, R. NatureWorks, LLC. Personal communication, 2008. (35) O’Connell, J. P.; Haile, J. M. Thermodynamics: Fundamentals for applications. Cambridge University Press: New York, 2005; p 82. (36) Boyd, R. H. Glass transition temperatures from molecular dynamics simulations. Trends Polym. Sci. 1996, 4 (1), 12–17. (37) Han, J.; Gee, R. H.; Boyd, R. H. Glass-transition temperatures of polymers from molecular-dynamics simulations. Macromolecules 1994, 27 (26), 7781–7784. (38) Rigby, D.; Roe, R. J. Molecular-dynamics simulation of polymer liquid and glass. I. Glass-transition. J. Chem. Phys. 1987, 87 (12), 7285– 7292. (39) Buchholz, J.; Paul, W.; Varnik, F.; Binder, K. Cooling rate dependence of the glass transition temperature of polymer melts: Molecular dynamics study. J. Chem. Phys. 2002, 117 (15), 7364–7372. (40) Soldera, A.; Metatla, N. Glass transition of polymers: Atomistic simulation versus experiments. Phys. Rev. E: Stat. Nonlin. Soft Matter Phys. 2006, 74 (61), 061803/1–061803/6. (41) Sperling, L. H. Introduction to physical polymer science, 4th ed.; Wiley: Hoboken, NJ, 2006. (42) 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 (14), 3701–3707. (43) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the opls-aa force field. J. Comput. Chem. 2001, 22 (13), 1340–1352. (44) Zhang, J.; Liang, Y.; Yan, J. Z.; Lou, J. Z. Study of the molecular weight dependence of glass transition temperature for amorphous poly(l-lactide) by molecular dynamics simulation. Polymer 2007, 48 (16), 4900–4905.

3767


Development of Force Field Parameters for Molecular Simulation of

Recommend Documents