- Email: [email protected]
Solubility of water in polymers—atomistic simulations
Computational and Theoretical Polymer Science 11 (2001) 49±55
www.elsevier.nl/locate/ctps
Solubility of water in polymersÐatomistic simulations B. Nick a, U.W. Suter b,* a
German Wool Research Institute at the Aachen University of Technology, Veltmanplatz 8, D-52062 Aachen, Germany b Department of Materials, Institute of Polymers, ETH, CH-8092 ZuÈrich, Switzerland Received 1 December 1998; received in revised form 16 June 1999; accepted 27 June 1999
Abstract The combination of the thermodynamic-integration approach and Widom's particle-insertion method is shown to be appropriate to determine the excess chemical potential of water in dense, amorphous polymer microstructures from MD simulation at an atomistically detailed level. The two-step method is applied to bisphenol-A±polycarbonate (BPA±PC) and polyvinylalcohol (PVA). The results are compared to previously published calculations on water sorption of various polyamides and show the applicability of the two-step method for the calculation of the excess chemical potential of water in a variety of polymer materials to obtain an estimate of their water sorption behavior. q 2000 Elsevier Science Ltd. All rights reserved. Keywords: Sorption; Widom's particle-insertion method; Polyamides
1. Introduction Sorption of small molecules by polymeric materials and the solubility of polymers in low-molecular-weight compounds are topics of long-standing interest. Especially the sorption of water in polymer systems is of great practical importance for nearly every application of these materials as evidenced, e.g. by the remarkable dependence of the mechanical properties of some polymers on water sorption. Further, understanding the sorption behavior is important in many areas of technology such as the design of barrier and separation membranes, packaging, and polymer processing. Molecular modeling provides a promising tool for investigation of the interplay between small molecules and the polymer microstructure and, in principle, opens up the possibility to calculate diffusion and solubility parameters. However, a detailed analysis of the sorption behavior requires atomistically detailed structures of the polymer. Due to limitations in computer power and poor methodological approaches, calculations of diffusion coef®cients and solubility parameters were limited to small, non-polar gases or to structural investigations of cavity sizes and distributions that are commonly thought to in¯uence the solubility of small molecules [1±19]. These calculations correctly describe the relative ordering of solubilities among different gases, but the quantitative values exceed the corresponding * Corresponding author. Fax: 1 41-1-632-1096. E-mail address: [email protected] (U.W. Suter)
experimental values by 1±2 orders of magnitude (except for the simplest compounds) [4,12]. Despite the fact that more sophisticated methods for polymer structure generation seem to improve the results [20], the extension of solubility calculations to larger molecules or molecules exhibiting long-range interactions has not yet succeeded satisfactorily. However, recent simulations on polyamide/water systems have shown that it is possible to calculate the equilibrium water content of different polyamides [21,22]. In these simulations, a new two-step method for the calculation of the excess chemical potential of water in dense polymer microstructures is introduced; this method successfully combines the thermodynamic integration approach and Widom's particle-insertion method to calculate the excess chemical potential of water in dense polymer microstructures. After comparison of the excess chemical potential of the dissolved-water and pure-water, an estimation of the equilibrium water content of the polymer material is obtained. A detailed description of the method is given in Ref. [21]. It is the aim of this study to demonstrate the general applicability of the two-step method. In Ref. [22] it was shown that calculation of the excess chemical potential of water in dense, amorphous polyamide microstructures is possible and gives an estimation of the equilibrium water content. Here, we extend our study to other types of polymer materials. The calculations are performed for bisphenol-A± polycarbonate (BPA±PC) and polyvinylalcohol (PVA), two polymers with completely different water sorption behavior.
1089-3156/01/$ - see front matter q 2000 Elsevier Science Ltd. All rights reserved. PII: S 1089-315 6(99)00061-6
50
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Table 1 The polymer/water model systems investigated, T 300 K
BPA±PC, X 30 rinitial 0:5 g cm23
PVA, X 100 rinitial 1:27 g cm23
Water content (wt%)
No. of water molecules
Ê ) (before equi.) Box edge (A
Cutoff correction pressure (bar)
2 5 10 15 20 30 10 20 30 50
9 22 47 75 106 181 27 61 105 245
29.58 29.87 30.42 31.01 31.64 33.07 18.56 19.31 20.20 22.59
1603 1603 1603 1602 1602 1602 2116 1757 1762 1771
The results are compared to those from previous studies carried out for polyamide(6) (PA-6) and polyamide(12) (PA-12) [22] to demonstrate the broad practicability of this method for water sorption estimations. 2. Simulation details The calculations on various BPA±PC/water and PVA/ water systems were performed on IBM RS 6000 or SGI indigo workstations, using a standard molecular modeling package (Insight 3.0.0/Discover 2.9.7, MSI, 1996) [23]. The Verlet leapfrog method [24] was used to integrate the equations of motions with a time-step of 1 fs. All MD simulations were carried out at 300 K using Berendsen's thermostat [25] to control temperature during data sampling. The pcff force ®eld was used to describe the atom±atom interactions at fully atomistic level, including all internal degrees-of-freedom for the water atoms as well as the polymer atoms. The non-bonded interactions, the Lennard± Jones and Coulombic interactions, were truncated at Ê Ða typical cut off distance for ªamorphous cellº 8.5 A Ê. simulationsÐusing a ®fth-order spline from 7.5 to 8.5 A Periodic continuation conditions and the minimum image convention applied to all simulations. The initial polymer/water systems were generated with the ªamorphous cellº procedure [23]; an RIS-based scheme determines the polymer backbone structure [26] and the water molecules are inserted randomly. The initial densities of the PVA/water systems were set at the experimentally determined density of the dry amorphous polymer rPVA;dry 1:27 g cm23 [28]. For the BPA±PC/water systems it is not possible to use the ªamorphous cellº procedure starting form the experimental density rBPA±PC;dry 1:2 g cm23 [28]; the model-generation algorithm stops most of the time because of ring catenation, especially for low-water content. To overcome this problem, the initial densities of the BPAPC/water systems were set to r 0:5 g cm23 and these structures were then densi®ed during the beginning of the equilibration. All model boxes were built of one polymer chain and a speci®ed number of water molecules according to the desired water content of the
systems. The number of monomeric units of the polymer chains (the length was chosen to maintain the ®nal box edge ÊÐ length of all systems investigated at smaller than 25 A resulting in reasonable computing times on current workstations) and the number of water molecules of all systems are given in Table 1. Because of the larger statistical error of the systems containing fewer water molecules, three independent structures each were generated of the BPAPC system containing 2 wt% water. All initial polymer/water structures were ®rst equilibrated to adapt the densities to the in¯uence of the water molecules. To this end, the potential energy with respect to the coordinates of all system atoms was ®rst subjected to 1000 steps of steepest-descent minimization. The PVA structures were then MD-equilibrated for 50 ps under constant-NVT conditions to pre-relax the polymer con®guration. Further equilibration and density adjustment were then achieved by an additional 550 ps constant-NpT MD simulation using Berendsen's method [25], which couples the system to a pressure ªbathº to maintain the pressure at a certain target. An external pressure was applied during the constant-NpT runs to compensate for the longe-range Lennard-Jones forces neglected by the short cutoff [27]. The BPA±PC systems were equilibrated for 500 ps with constant-NpT MD simulation. A fast compression of the low starting density was achieved by the external pressure calculated according to a system density of r 1:2 g cm23 (the density of dry BPA±PC). The correction pressure values are also given in Table 1. The last con®guration of each equilibration run was utilized as the initial con®guration for all MD runs of the thermodynamic integration scheme. This means that data sampling at all intermediate points of the thermodynamic integration (applying the l -dependent force-®elds) started from identical con®gurations. To adapt the con®gurations to these force-®elds, the structures were again equilibrated for 50 ps by NVT MD; afterwards, 150 (in case of BPA±PC) and 80 (in case of PVA) con®gurations of the dynamic trajectories were sampled at 1 ps intervals, employing the l -dependent force-®elds. These con®gurations are used to carry out thermodynamic integration and to calculate the difference in Helmholtz free energy due to
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
51
3. Results and discussion 3.1. Density
Fig. 1. Density vs time during equilibration of a BPA±PC and a PVA system containing 20 wt% water each. The curves show instantaneous values sampled at 1 ps intervals during NpT-MD simulation.
the non-bonded-interaction scaling of the water molecules. A detailed overview of the entire simulation scheme applied to every polymer/water structure is given in Ref. [22]. The l -dependent force parameters of the water atoms that are valid during the MD simulations of the thermodynamic integration can be found in Ref. [21]. The con®gurations collected during the sampling run with the l 0 force-®eldÐat this stage the water is called ªnoble waterº considering its spherical shape and small nonbonded interactions (although the internal degrees-of-freedom of real water are retained)Ðthe polymer/noble water systems, are further utilized to calculate the excess chemical potential of noble water in the amorphous polymer microstructure by means of Widom's particle-insertion method. As already described in our ®rst paper concerning the two step-method [21], we generate a 30 £ 30 £ 30 grid inside the model boxes and perform a noble-water insertion on all grid points to obtain the excess chemical potential of the noble water dissolved in the dense polymer matrix.
Equilibration and density adjustment of the initial polymer/water systems is achieved through MD simulations under constant NpT-conditions. The equilibration progress of two systems is monitored in Fig. 1. Because the BPA± PC/water are generated with a low density r 0:5 g cm23 ; the system density in Fig. 1 shows a strong increase during the ®rst 50 ps NpT-MD simulation, before it settles down towards a nearly constant level up to the end of the equilibration after 500 ps. The PVA/water systems are generated with the density of dry amorphous PVA. The equilibration is started with 50 ps under constant NVT-conditions (causing the constant-density stretch at the beginning of the PVAcurve in Fig. 1) to avoid strong simulation-box deformations. After this initial time, the density decreases abruptly when changing to NpT-conditions, accounting for the water in¯uence, and settles down to a constant level for the rest of the equilibration. Although our equilibration process lasts for 500 and 600 ps, respectively, it is worth emphasizing that the equilibrated polymer/water systems are still near their initial con®guration. After this equilibration time the systems seem to have reached a local minima in their con®guration space; the system densities do not show any drifts larger than 3 £ 1024 g cm21 ps21 during the last 100 ps of the equilibration (this criterion was used to follow the equilibration progress). Even after much longer equilibration time no signi®cant changes can be expected. However, the calculated properties, i.e. the excess chemical potential of water in the polymer microstructures, might depend on the speci®c polymer con®guration. This has to be taken into account when interpreting the results and will be investigated for the BPA±PC/water system containing 2 wt% water. The water in¯uence on the density after equilibration is shown in Fig. 2 for all systems and compared to the density that is expected when applying a simple linear combination rule (by weight fraction) for the estimation of the density of mixed systems. Even though no de®nitive correlation can be established with so few data points, the linear combination rule seems to apply for all systems investigated within the precision characteristic for density data calculated by MD simulations. However, for up to 30 wt% water in the BPA± PC microstructure, no signi®cant decrease of the system density during the equilibration time can be observed. The individual system densities after equilibration are given in Table 2 and are used for the subsequent Widom's particleinsertion. 3.2. Water sorption/solubility
Fig. 2. Density after equilibration as a function of the water content of the polymer/water systems investigated (lines are for a linear combination rule, by weight fraction of water). The dry-polymer densities were taken from the literature rBPA±PC 1:2 g cm23 ; and rPVA 1:27 g cm23 [28].
Estimation of the water sorption or solubility of polymer materials is possible by comparison of the excess chemical potential of water dissolved in the polymer matrix with that
52
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Table 2 Results for the polymer/water structures, T 300 K
BPA±PC
PVA
H2 O
Water content (wt%)
DATI =NW (kJ mol 23)
r (g cm 23)
mWI (kJ mol 23)
m ex (kJ mol 23)
2 2 2 2 5 10 15 20 30 10 20 30 50 100
218.3 220.5 218.6
1.18 1.18 1.19 1.18 1.18 1.20 1.18 1.18 1.15 1.25 1.20 1.20 1.16 1.00
20.9 20.1 20.1
219.2 220.6 218.7 219.5 219.8 219.4 221.0 220.7 222.0 224.3 222.7 224.1 224.6 223.8
219.7 219.9 220.7 220.6 221.4 228.7 225.5 226.2 224.5 222.2
of pure water. Calculation of the excess chemical potential becomes feasible by applying our two-step method that combines thermodynamic integration and Widom's particle-insertion as introduced in Ref. [21]. An overview of the calculated results for the excess chemical potential of water in the BPA±PC and PVA systems investigated is given in Table 2. The main contribution to the excess chemical potential seems to originate from the scaling of the non-bonded interactionsÐespecially the Coulombic interactions that dominate in polymer/water systemsÐduring thermodynamic integration. The Coulombic interactions of the water molecules are completely eliminated during the thermodynamic integration and cause a large change in Helmholtz energy. The remaining non-bonded interactions of the water molecules represent the Lennard-Jones interactions of the wateroxygen atoms. The contribution of these interactions to the excess chemical potential is calculated by Widom's particle-insertion and gives only a minor contribution to the Helmholtz free energy of the water molecules. Although the values of the intermediate results calculated
Fig. 3. Difference in Helmholtz energy between water and noble water dissolved in different polymers, calculated with thermodynamic integration.
20.1 10.5 20.3 20.1 20.6 14.4 12.8 12.1 20.1 21.6(0.1)
for BPA±PC exhibit a general similarity in size and ratio to the values calculated in our previous study [22] for polyamide/water systems, they do not evolve equally. The difference in Helmholtz free energy between water and noble water calculated by thermodynamic integration slightly decreases with increasing water content for the BPA±PC/ water systems (see Fig. 3), even though the effect is not very pronounced and the BPA±PC values do not differ more than kT from the correspondent PA-12 values. The functional form of 2Upot l=2l does not change when compared to the polyamide/water or pure water systems of Ref. [21]. As demonstrated in Fig. 4, 2Upot l=2l still decreases nearly linearly with increasing l , con®rming that the Coulombic interactions dominate the potential energy also in BPA±PC/ water systems for higher l -values. In BPA±PC we do not ®nd the strong systematic decrease of the excess chemical potential of dissolved noble water calculated with Widom's particle-insertion observed for all polyamide/water systems investigated earlier (Fig. 5). The excess chemical potential of noble water in systems of identical composition shows larger variations than changes caused by increasing water content. However, if the ®rst BPA±PC/system containing 2 wt% water of Table 2 is
Fig. 4. Development of the derivative of Upot l with respect to the coupling parameter l for different BPA±PC/water systems.
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Fig. 5. Excess chemical potential of noble water in various polymer/noble water systems, calculated by Widom's particle-insertion method, plotted against the number density of the water molecules. Lines are linear ®ts used as guides for the eye.
regarded to be an outlier, a slightly decreasing regression line is also obtained for the BPA±PC structures; the stronger in¯uence of the speci®c system rather than the increasing water content (causing more `empty' space in the system) was already pointed out in Ref. [22]. In contrast to the BPA±PC/water systems, water dissolved in PVA shows a completely different behavior. Both the intermediate results evolve with increasing water content as for PA-6 [22]. The Helmholtz energy difference between water and noble water dissolved in the PVA microstructure increases with increasing water content (Fig. 3); the excess chemical potential of noble water decreases with increasing water content (Fig. 5). Therefore the development of both the intermediate results corresponds directly to the behavior already observed for PA-6. The sorption behavior of the polymer systems submerged in neat water can be determined by looking at the difference between the excess chemical potential of water dissolved in
53
the polymer microstructures and that of pure water. The excess chemical potential of water dissolved in the polymer structures is calculated as the sum of the intermediate results (thermodynamic integration and Widom's particle-insertion method) listed in Table 2. The excess chemical potential of pure water, the reference state, was taken from Ref. [21] as 21 mex w 223:8 kJ mol ; calculated using exactly the same method and assumptions as in the current study. The calculated values for the difference in excess chemical potential between water dissolved in various polymer microstructures and pure water are shown in Fig. 6. Fig. 6 also replicates some values calculated for PA-6 and PA-12 from Ref. [22] for comparison. For all BPA±PC/water systems investigated, the difference between the excess chemical potential of water dissolved in the polymer microstructure and pure water is positive, indicating that the BPA±PC systems investigated will absorb no water from a wet environment. The values lie in the same range that was already calculated for PA12/ water systems and lead to the analogous conclusion: the smallest water content investigated (2 wt% water) is already higher than the equilibrium water content of BPA±PC. This result derived for BPA±PC is in reasonable agreement with experimental values that measure the saturation point of water in BPA±PC at 0.3 wt% [28]. This low-water content cannot be simulated with the box-size restrictions given by the available computer time; our two-step method requires at least a few water molecules inside the simulation box to give statistically acceptable results. It is even not clear that the experimental water content is mainly determined by the sorption capability of the polymer material and does not re¯ect water sorption in structural irregularities of the specimen (some kind of Langmuir sorption). In contrast, the sorption behavior of PVA turns out to be nearly independent of the water concentration of the microstructure in our simulations. The ®nal values for the excess chemical potential of water dissolved in the four PVA/water systems investigated do not signi®cantly differ from that of pure water; in all cases, the difference is less than 0.5 kT. From our calculations, we expect amorphous PVA to be soluble in water, in agreement with experiments. We can also deduce the concentration dependence of the polymer/solvent interactions from our simulations. Polymer/solvent interactions have been described successfully by the Flory±Huggins theory [29,30]. This theory allows the estimation of the excess chemical potential of a solvent dissolved in a polymer matrix from Dmex RTln F1 1 12F1 1x 12F1 2
Fig. 6. Calculated difference between the excess chemical potential of water dissolved in amorphous BPA±PC and PVA microstructures and that of pure water at 300 K. Polyamide/water results are added from Ref. [22] for comparison. The curves shown are simply guides for the eye.
with F 1 being the volume fraction of the solvent (water) in the polymer matrix, x the parameter related to the interaction Helmholtz energy of the ®rst neighbors, R the universal gas constant and T the temperature. Rearrangement of this equation and the utilization of our simulation results allows the determination of the interaction parameter x , taking into
54
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Fig. 7. Polymer/solvent interaction parameter x , calculated from the simulated differential excess chemical potentials, as a function of the water content of the systems. Lines are linear ®ts to the data points.
consideration the possibility of a functional dependence of x on the solvent concentration. For this, the calculated interaction parameter x is plotted vs. the water content of the polymer microstructures in Fig. 7. For both polymers (BPA±PC and PVA), our calculations show a concentration-dependent interaction parameter. If one used a simple linear ®t for the description of x F1 ; one would arrive at xBPA±PC F1 4:726:6F1 and xPVA F1 2:625:9F1 ; respectively. Using these x F1 parameters, one is able to draw smooth curves consistent with the Flory±Huggins theory into Fig. 6 to calculate the sorption limits for the speci®c polymers investigated. For BPA±PC and PVA this does not make much sense because of the very small number that would be calculated in case of BPA±PCÐwhich is surely outside the precision of our simulationsÐand the strong in¯uence of scatter of the few simulation points in the case of PVA (only 80 frames were sampled for every l dependent force®eld). Nevertheless, we want to point out the possibility of investigation of the concentration dependence of the polymer/solvent interaction parameter with our calculations. 4. Conclusion The calculations presented in this paper for BPA±PC and PVA and the earlier results for different polyamide systems [22] demonstrate clearly that our new two-step methodÐ the combination of thermodynamic integration and the Widom's particle-insertion methodÐgives a good estimate of the excess chemical potential of water in dense polymer/ microstructures. From our calculations we can obtain an approximation of the degree of sorption of water into the amorphous microstructure. This estimation is not restricted to speci®c polymeric materials. The quality of such calculations only depends on the energy model used, i.e. on the force-®eld for the atom±atom interactions during the MD simulations. The increasing reliability of the force-®elds for
polymeric materials in recent years makes the approach presented highly promising for future materials characterization. As already mentioned, the differential excess chemical potentials are small (within a few kJ mol 21); nevertheless, we want to stress that we can reproducibly distinguish between different sorption behaviors for different polymer materials. Conceding that calculations of the equilibrium water content are restricted by the number of water molecules in the model systemsÐonly a few water molecules inside the system lead to less expedient statistics during data samplingÐour method allows the determination of equilibrium sorption data of polymers at experimental densities in a thermodynamically sound way. Besides, our two-step method allows us to obtain structural information from the MD runs, making it superior to some other methods that produce biased trajectories, e.g. umbrella sampling. The general application of this two-step method for the calculation of excess chemical potentials of small guests in polymer microstructures seems to be straightforward. A prerequisite is that an appropriate scaling function of the non-bonded interaction is found so that Widom's particleinsertion method becomes feasible. In theory, after addition of the two intermediate results, the in¯uence of the scaling function should vanish. However, it turns out that ®nding an ef®cient scaling law can be dif®cult; at present, we cannot give a general recipe for the construction of scaling functions. Also, sampling must be ef®cient (an inherent problem of MD simulations) and this might limit the applicability of the method; nevertheless, the approach presented is promising and the assistance of more sophisticated sampling methods [31] might help to extend the investigations to more complex compounds. Acknowledgements We gratefully acknowledge support for B.N. from the Deutsche Forschungs-gemeinschaft under Grant no. DFG Kn 417/12 and general ®nancial support by EMS-Chemie AG, Domat/Ems, Switzerland. Generous amounts of computer time was provided by the Competence Centre for Computational Chemistry (C4) of the ETH ZuÈrich. We also acknowledge valuable discussions with Andrei A. Gusev. References [1] [2] [3] [4] [5] [6]
Shing KS, Gubbins KE. Mol Phys 1981;43:717. Takeuchi H, Okazaki K. J Chem Phys 1990;92:5643. Takeuchi H, Roe R-J, Mark JE. J Chem Phys 1990;93:9042. Mueller-Plathe F. Macromolecules 1991;24:6475. Gusev AA, Suter UW. Phys Rev A 1991;43:6488. Sok RM, Berendsen HJC, van Gunsteren WF. J Chem Phys 1992;96:4699. [7] de Pablo JJ, Laso M, Suter UW. J Chem Phys 1992;96:6157. [8] de Pablo JJ, Laso M, Suter UW. Macromolecules 1993;26:6180.
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55 [9] Pant PV, Boyd RH. J Chem Phys 1993;98:9895. [10] Gusev AA, Arizzi S, Suter UW, Moll DJ. J Chem Phys 1993;99: 2221. [11] Gusev AA, Suter UW. J Chem Phys 1993;99:2228. [12] Mueller-Plathe F, Rogers SC, van Gunsteren WF. J Chem Phys 1993;98:9895. [13] Han J, Boyd RH. Macromolecules 1994;27:5365. [14] Tamai Y, Tanaka H, Nakanishi K. Macromolecules 1994;27:4498. [15] Tamai Y, Tanaka H, Nakanishi K. Macromolecules 1995;28:2544. [16] Gee RH, Boyd RH. Polymer 1995;36:1435. [17] Gusev AA, Suter UW, Moll DJ. Macromolecules 1995;28:2582. [18] Ralston ARK, Denton DD, Bicerano J, Moll DJ. Comput Polym Sci 1996;6:15. [19] Bicerano J, Moll DJ. Comput Polym Sci 1996;6:117. [20] van der Vegt NFA, Briels WJ, Wessling M, Strathmann H. J Chem Phys 1997;105:8849.
55
[21] Knopp B, Suter UW, Gusev AA. Macromolecules 1997;30:6107± 6113. [22] Knopp B, Suter UW. Macromolecules 1997;30:6114±6119. [23] Insight 3.0.0. Discover 2.9.7. Amorphous Cell 8.0.0, MSI, San Diego, CA, 1996. [24] Verlet L. Phys Rev 1967;159:98. [25] Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. J Chem Phys 1984;81:3684. [26] Theodorou DN, Suter UW. Macromolecules 1985;18:1467. [27] Allen MP, Tildesley DJ. Computer simulation of liquids, Oxford University Press, 1987. [28] Brandrup J, Immergut EH, editors. Polymer handbook 3. New York: Wiley, 1989. [29] Flory PJ. J Chem Phys 1942;10:51. [30] Huggins ML. J Chem Phys 1941;9:440. [31] Leontidis E, Suter UW. Mol Phys 1994;83:489.
www.elsevier.nl/locate/ctps
Solubility of water in polymersÐatomistic simulations B. Nick a, U.W. Suter b,* a
German Wool Research Institute at the Aachen University of Technology, Veltmanplatz 8, D-52062 Aachen, Germany b Department of Materials, Institute of Polymers, ETH, CH-8092 ZuÈrich, Switzerland Received 1 December 1998; received in revised form 16 June 1999; accepted 27 June 1999
Abstract The combination of the thermodynamic-integration approach and Widom's particle-insertion method is shown to be appropriate to determine the excess chemical potential of water in dense, amorphous polymer microstructures from MD simulation at an atomistically detailed level. The two-step method is applied to bisphenol-A±polycarbonate (BPA±PC) and polyvinylalcohol (PVA). The results are compared to previously published calculations on water sorption of various polyamides and show the applicability of the two-step method for the calculation of the excess chemical potential of water in a variety of polymer materials to obtain an estimate of their water sorption behavior. q 2000 Elsevier Science Ltd. All rights reserved. Keywords: Sorption; Widom's particle-insertion method; Polyamides
1. Introduction Sorption of small molecules by polymeric materials and the solubility of polymers in low-molecular-weight compounds are topics of long-standing interest. Especially the sorption of water in polymer systems is of great practical importance for nearly every application of these materials as evidenced, e.g. by the remarkable dependence of the mechanical properties of some polymers on water sorption. Further, understanding the sorption behavior is important in many areas of technology such as the design of barrier and separation membranes, packaging, and polymer processing. Molecular modeling provides a promising tool for investigation of the interplay between small molecules and the polymer microstructure and, in principle, opens up the possibility to calculate diffusion and solubility parameters. However, a detailed analysis of the sorption behavior requires atomistically detailed structures of the polymer. Due to limitations in computer power and poor methodological approaches, calculations of diffusion coef®cients and solubility parameters were limited to small, non-polar gases or to structural investigations of cavity sizes and distributions that are commonly thought to in¯uence the solubility of small molecules [1±19]. These calculations correctly describe the relative ordering of solubilities among different gases, but the quantitative values exceed the corresponding * Corresponding author. Fax: 1 41-1-632-1096. E-mail address: [email protected] (U.W. Suter)
experimental values by 1±2 orders of magnitude (except for the simplest compounds) [4,12]. Despite the fact that more sophisticated methods for polymer structure generation seem to improve the results [20], the extension of solubility calculations to larger molecules or molecules exhibiting long-range interactions has not yet succeeded satisfactorily. However, recent simulations on polyamide/water systems have shown that it is possible to calculate the equilibrium water content of different polyamides [21,22]. In these simulations, a new two-step method for the calculation of the excess chemical potential of water in dense polymer microstructures is introduced; this method successfully combines the thermodynamic integration approach and Widom's particle-insertion method to calculate the excess chemical potential of water in dense polymer microstructures. After comparison of the excess chemical potential of the dissolved-water and pure-water, an estimation of the equilibrium water content of the polymer material is obtained. A detailed description of the method is given in Ref. [21]. It is the aim of this study to demonstrate the general applicability of the two-step method. In Ref. [22] it was shown that calculation of the excess chemical potential of water in dense, amorphous polyamide microstructures is possible and gives an estimation of the equilibrium water content. Here, we extend our study to other types of polymer materials. The calculations are performed for bisphenol-A± polycarbonate (BPA±PC) and polyvinylalcohol (PVA), two polymers with completely different water sorption behavior.
1089-3156/01/$ - see front matter q 2000 Elsevier Science Ltd. All rights reserved. PII: S 1089-315 6(99)00061-6
50
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Table 1 The polymer/water model systems investigated, T 300 K
BPA±PC, X 30 rinitial 0:5 g cm23
PVA, X 100 rinitial 1:27 g cm23
Water content (wt%)
No. of water molecules
Ê ) (before equi.) Box edge (A
Cutoff correction pressure (bar)
2 5 10 15 20 30 10 20 30 50
9 22 47 75 106 181 27 61 105 245
29.58 29.87 30.42 31.01 31.64 33.07 18.56 19.31 20.20 22.59
1603 1603 1603 1602 1602 1602 2116 1757 1762 1771
The results are compared to those from previous studies carried out for polyamide(6) (PA-6) and polyamide(12) (PA-12) [22] to demonstrate the broad practicability of this method for water sorption estimations. 2. Simulation details The calculations on various BPA±PC/water and PVA/ water systems were performed on IBM RS 6000 or SGI indigo workstations, using a standard molecular modeling package (Insight 3.0.0/Discover 2.9.7, MSI, 1996) [23]. The Verlet leapfrog method [24] was used to integrate the equations of motions with a time-step of 1 fs. All MD simulations were carried out at 300 K using Berendsen's thermostat [25] to control temperature during data sampling. The pcff force ®eld was used to describe the atom±atom interactions at fully atomistic level, including all internal degrees-of-freedom for the water atoms as well as the polymer atoms. The non-bonded interactions, the Lennard± Jones and Coulombic interactions, were truncated at Ê Ða typical cut off distance for ªamorphous cellº 8.5 A Ê. simulationsÐusing a ®fth-order spline from 7.5 to 8.5 A Periodic continuation conditions and the minimum image convention applied to all simulations. The initial polymer/water systems were generated with the ªamorphous cellº procedure [23]; an RIS-based scheme determines the polymer backbone structure [26] and the water molecules are inserted randomly. The initial densities of the PVA/water systems were set at the experimentally determined density of the dry amorphous polymer rPVA;dry 1:27 g cm23 [28]. For the BPA±PC/water systems it is not possible to use the ªamorphous cellº procedure starting form the experimental density rBPA±PC;dry 1:2 g cm23 [28]; the model-generation algorithm stops most of the time because of ring catenation, especially for low-water content. To overcome this problem, the initial densities of the BPAPC/water systems were set to r 0:5 g cm23 and these structures were then densi®ed during the beginning of the equilibration. All model boxes were built of one polymer chain and a speci®ed number of water molecules according to the desired water content of the
systems. The number of monomeric units of the polymer chains (the length was chosen to maintain the ®nal box edge ÊÐ length of all systems investigated at smaller than 25 A resulting in reasonable computing times on current workstations) and the number of water molecules of all systems are given in Table 1. Because of the larger statistical error of the systems containing fewer water molecules, three independent structures each were generated of the BPAPC system containing 2 wt% water. All initial polymer/water structures were ®rst equilibrated to adapt the densities to the in¯uence of the water molecules. To this end, the potential energy with respect to the coordinates of all system atoms was ®rst subjected to 1000 steps of steepest-descent minimization. The PVA structures were then MD-equilibrated for 50 ps under constant-NVT conditions to pre-relax the polymer con®guration. Further equilibration and density adjustment were then achieved by an additional 550 ps constant-NpT MD simulation using Berendsen's method [25], which couples the system to a pressure ªbathº to maintain the pressure at a certain target. An external pressure was applied during the constant-NpT runs to compensate for the longe-range Lennard-Jones forces neglected by the short cutoff [27]. The BPA±PC systems were equilibrated for 500 ps with constant-NpT MD simulation. A fast compression of the low starting density was achieved by the external pressure calculated according to a system density of r 1:2 g cm23 (the density of dry BPA±PC). The correction pressure values are also given in Table 1. The last con®guration of each equilibration run was utilized as the initial con®guration for all MD runs of the thermodynamic integration scheme. This means that data sampling at all intermediate points of the thermodynamic integration (applying the l -dependent force-®elds) started from identical con®gurations. To adapt the con®gurations to these force-®elds, the structures were again equilibrated for 50 ps by NVT MD; afterwards, 150 (in case of BPA±PC) and 80 (in case of PVA) con®gurations of the dynamic trajectories were sampled at 1 ps intervals, employing the l -dependent force-®elds. These con®gurations are used to carry out thermodynamic integration and to calculate the difference in Helmholtz free energy due to
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
51
3. Results and discussion 3.1. Density
Fig. 1. Density vs time during equilibration of a BPA±PC and a PVA system containing 20 wt% water each. The curves show instantaneous values sampled at 1 ps intervals during NpT-MD simulation.
the non-bonded-interaction scaling of the water molecules. A detailed overview of the entire simulation scheme applied to every polymer/water structure is given in Ref. [22]. The l -dependent force parameters of the water atoms that are valid during the MD simulations of the thermodynamic integration can be found in Ref. [21]. The con®gurations collected during the sampling run with the l 0 force-®eldÐat this stage the water is called ªnoble waterº considering its spherical shape and small nonbonded interactions (although the internal degrees-of-freedom of real water are retained)Ðthe polymer/noble water systems, are further utilized to calculate the excess chemical potential of noble water in the amorphous polymer microstructure by means of Widom's particle-insertion method. As already described in our ®rst paper concerning the two step-method [21], we generate a 30 £ 30 £ 30 grid inside the model boxes and perform a noble-water insertion on all grid points to obtain the excess chemical potential of the noble water dissolved in the dense polymer matrix.
Equilibration and density adjustment of the initial polymer/water systems is achieved through MD simulations under constant NpT-conditions. The equilibration progress of two systems is monitored in Fig. 1. Because the BPA± PC/water are generated with a low density r 0:5 g cm23 ; the system density in Fig. 1 shows a strong increase during the ®rst 50 ps NpT-MD simulation, before it settles down towards a nearly constant level up to the end of the equilibration after 500 ps. The PVA/water systems are generated with the density of dry amorphous PVA. The equilibration is started with 50 ps under constant NVT-conditions (causing the constant-density stretch at the beginning of the PVAcurve in Fig. 1) to avoid strong simulation-box deformations. After this initial time, the density decreases abruptly when changing to NpT-conditions, accounting for the water in¯uence, and settles down to a constant level for the rest of the equilibration. Although our equilibration process lasts for 500 and 600 ps, respectively, it is worth emphasizing that the equilibrated polymer/water systems are still near their initial con®guration. After this equilibration time the systems seem to have reached a local minima in their con®guration space; the system densities do not show any drifts larger than 3 £ 1024 g cm21 ps21 during the last 100 ps of the equilibration (this criterion was used to follow the equilibration progress). Even after much longer equilibration time no signi®cant changes can be expected. However, the calculated properties, i.e. the excess chemical potential of water in the polymer microstructures, might depend on the speci®c polymer con®guration. This has to be taken into account when interpreting the results and will be investigated for the BPA±PC/water system containing 2 wt% water. The water in¯uence on the density after equilibration is shown in Fig. 2 for all systems and compared to the density that is expected when applying a simple linear combination rule (by weight fraction) for the estimation of the density of mixed systems. Even though no de®nitive correlation can be established with so few data points, the linear combination rule seems to apply for all systems investigated within the precision characteristic for density data calculated by MD simulations. However, for up to 30 wt% water in the BPA± PC microstructure, no signi®cant decrease of the system density during the equilibration time can be observed. The individual system densities after equilibration are given in Table 2 and are used for the subsequent Widom's particleinsertion. 3.2. Water sorption/solubility
Fig. 2. Density after equilibration as a function of the water content of the polymer/water systems investigated (lines are for a linear combination rule, by weight fraction of water). The dry-polymer densities were taken from the literature rBPA±PC 1:2 g cm23 ; and rPVA 1:27 g cm23 [28].
Estimation of the water sorption or solubility of polymer materials is possible by comparison of the excess chemical potential of water dissolved in the polymer matrix with that
52
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Table 2 Results for the polymer/water structures, T 300 K
BPA±PC
PVA
H2 O
Water content (wt%)
DATI =NW (kJ mol 23)
r (g cm 23)
mWI (kJ mol 23)
m ex (kJ mol 23)
2 2 2 2 5 10 15 20 30 10 20 30 50 100
218.3 220.5 218.6
1.18 1.18 1.19 1.18 1.18 1.20 1.18 1.18 1.15 1.25 1.20 1.20 1.16 1.00
20.9 20.1 20.1
219.2 220.6 218.7 219.5 219.8 219.4 221.0 220.7 222.0 224.3 222.7 224.1 224.6 223.8
219.7 219.9 220.7 220.6 221.4 228.7 225.5 226.2 224.5 222.2
of pure water. Calculation of the excess chemical potential becomes feasible by applying our two-step method that combines thermodynamic integration and Widom's particle-insertion as introduced in Ref. [21]. An overview of the calculated results for the excess chemical potential of water in the BPA±PC and PVA systems investigated is given in Table 2. The main contribution to the excess chemical potential seems to originate from the scaling of the non-bonded interactionsÐespecially the Coulombic interactions that dominate in polymer/water systemsÐduring thermodynamic integration. The Coulombic interactions of the water molecules are completely eliminated during the thermodynamic integration and cause a large change in Helmholtz energy. The remaining non-bonded interactions of the water molecules represent the Lennard-Jones interactions of the wateroxygen atoms. The contribution of these interactions to the excess chemical potential is calculated by Widom's particle-insertion and gives only a minor contribution to the Helmholtz free energy of the water molecules. Although the values of the intermediate results calculated
Fig. 3. Difference in Helmholtz energy between water and noble water dissolved in different polymers, calculated with thermodynamic integration.
20.1 10.5 20.3 20.1 20.6 14.4 12.8 12.1 20.1 21.6(0.1)
for BPA±PC exhibit a general similarity in size and ratio to the values calculated in our previous study [22] for polyamide/water systems, they do not evolve equally. The difference in Helmholtz free energy between water and noble water calculated by thermodynamic integration slightly decreases with increasing water content for the BPA±PC/ water systems (see Fig. 3), even though the effect is not very pronounced and the BPA±PC values do not differ more than kT from the correspondent PA-12 values. The functional form of 2Upot l=2l does not change when compared to the polyamide/water or pure water systems of Ref. [21]. As demonstrated in Fig. 4, 2Upot l=2l still decreases nearly linearly with increasing l , con®rming that the Coulombic interactions dominate the potential energy also in BPA±PC/ water systems for higher l -values. In BPA±PC we do not ®nd the strong systematic decrease of the excess chemical potential of dissolved noble water calculated with Widom's particle-insertion observed for all polyamide/water systems investigated earlier (Fig. 5). The excess chemical potential of noble water in systems of identical composition shows larger variations than changes caused by increasing water content. However, if the ®rst BPA±PC/system containing 2 wt% water of Table 2 is
Fig. 4. Development of the derivative of Upot l with respect to the coupling parameter l for different BPA±PC/water systems.
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Fig. 5. Excess chemical potential of noble water in various polymer/noble water systems, calculated by Widom's particle-insertion method, plotted against the number density of the water molecules. Lines are linear ®ts used as guides for the eye.
regarded to be an outlier, a slightly decreasing regression line is also obtained for the BPA±PC structures; the stronger in¯uence of the speci®c system rather than the increasing water content (causing more `empty' space in the system) was already pointed out in Ref. [22]. In contrast to the BPA±PC/water systems, water dissolved in PVA shows a completely different behavior. Both the intermediate results evolve with increasing water content as for PA-6 [22]. The Helmholtz energy difference between water and noble water dissolved in the PVA microstructure increases with increasing water content (Fig. 3); the excess chemical potential of noble water decreases with increasing water content (Fig. 5). Therefore the development of both the intermediate results corresponds directly to the behavior already observed for PA-6. The sorption behavior of the polymer systems submerged in neat water can be determined by looking at the difference between the excess chemical potential of water dissolved in
53
the polymer microstructures and that of pure water. The excess chemical potential of water dissolved in the polymer structures is calculated as the sum of the intermediate results (thermodynamic integration and Widom's particle-insertion method) listed in Table 2. The excess chemical potential of pure water, the reference state, was taken from Ref. [21] as 21 mex w 223:8 kJ mol ; calculated using exactly the same method and assumptions as in the current study. The calculated values for the difference in excess chemical potential between water dissolved in various polymer microstructures and pure water are shown in Fig. 6. Fig. 6 also replicates some values calculated for PA-6 and PA-12 from Ref. [22] for comparison. For all BPA±PC/water systems investigated, the difference between the excess chemical potential of water dissolved in the polymer microstructure and pure water is positive, indicating that the BPA±PC systems investigated will absorb no water from a wet environment. The values lie in the same range that was already calculated for PA12/ water systems and lead to the analogous conclusion: the smallest water content investigated (2 wt% water) is already higher than the equilibrium water content of BPA±PC. This result derived for BPA±PC is in reasonable agreement with experimental values that measure the saturation point of water in BPA±PC at 0.3 wt% [28]. This low-water content cannot be simulated with the box-size restrictions given by the available computer time; our two-step method requires at least a few water molecules inside the simulation box to give statistically acceptable results. It is even not clear that the experimental water content is mainly determined by the sorption capability of the polymer material and does not re¯ect water sorption in structural irregularities of the specimen (some kind of Langmuir sorption). In contrast, the sorption behavior of PVA turns out to be nearly independent of the water concentration of the microstructure in our simulations. The ®nal values for the excess chemical potential of water dissolved in the four PVA/water systems investigated do not signi®cantly differ from that of pure water; in all cases, the difference is less than 0.5 kT. From our calculations, we expect amorphous PVA to be soluble in water, in agreement with experiments. We can also deduce the concentration dependence of the polymer/solvent interactions from our simulations. Polymer/solvent interactions have been described successfully by the Flory±Huggins theory [29,30]. This theory allows the estimation of the excess chemical potential of a solvent dissolved in a polymer matrix from Dmex RTln F1 1 12F1 1x 12F1 2
Fig. 6. Calculated difference between the excess chemical potential of water dissolved in amorphous BPA±PC and PVA microstructures and that of pure water at 300 K. Polyamide/water results are added from Ref. [22] for comparison. The curves shown are simply guides for the eye.
with F 1 being the volume fraction of the solvent (water) in the polymer matrix, x the parameter related to the interaction Helmholtz energy of the ®rst neighbors, R the universal gas constant and T the temperature. Rearrangement of this equation and the utilization of our simulation results allows the determination of the interaction parameter x , taking into
54
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55
Fig. 7. Polymer/solvent interaction parameter x , calculated from the simulated differential excess chemical potentials, as a function of the water content of the systems. Lines are linear ®ts to the data points.
consideration the possibility of a functional dependence of x on the solvent concentration. For this, the calculated interaction parameter x is plotted vs. the water content of the polymer microstructures in Fig. 7. For both polymers (BPA±PC and PVA), our calculations show a concentration-dependent interaction parameter. If one used a simple linear ®t for the description of x F1 ; one would arrive at xBPA±PC F1 4:726:6F1 and xPVA F1 2:625:9F1 ; respectively. Using these x F1 parameters, one is able to draw smooth curves consistent with the Flory±Huggins theory into Fig. 6 to calculate the sorption limits for the speci®c polymers investigated. For BPA±PC and PVA this does not make much sense because of the very small number that would be calculated in case of BPA±PCÐwhich is surely outside the precision of our simulationsÐand the strong in¯uence of scatter of the few simulation points in the case of PVA (only 80 frames were sampled for every l dependent force®eld). Nevertheless, we want to point out the possibility of investigation of the concentration dependence of the polymer/solvent interaction parameter with our calculations. 4. Conclusion The calculations presented in this paper for BPA±PC and PVA and the earlier results for different polyamide systems [22] demonstrate clearly that our new two-step methodÐ the combination of thermodynamic integration and the Widom's particle-insertion methodÐgives a good estimate of the excess chemical potential of water in dense polymer/ microstructures. From our calculations we can obtain an approximation of the degree of sorption of water into the amorphous microstructure. This estimation is not restricted to speci®c polymeric materials. The quality of such calculations only depends on the energy model used, i.e. on the force-®eld for the atom±atom interactions during the MD simulations. The increasing reliability of the force-®elds for
polymeric materials in recent years makes the approach presented highly promising for future materials characterization. As already mentioned, the differential excess chemical potentials are small (within a few kJ mol 21); nevertheless, we want to stress that we can reproducibly distinguish between different sorption behaviors for different polymer materials. Conceding that calculations of the equilibrium water content are restricted by the number of water molecules in the model systemsÐonly a few water molecules inside the system lead to less expedient statistics during data samplingÐour method allows the determination of equilibrium sorption data of polymers at experimental densities in a thermodynamically sound way. Besides, our two-step method allows us to obtain structural information from the MD runs, making it superior to some other methods that produce biased trajectories, e.g. umbrella sampling. The general application of this two-step method for the calculation of excess chemical potentials of small guests in polymer microstructures seems to be straightforward. A prerequisite is that an appropriate scaling function of the non-bonded interaction is found so that Widom's particleinsertion method becomes feasible. In theory, after addition of the two intermediate results, the in¯uence of the scaling function should vanish. However, it turns out that ®nding an ef®cient scaling law can be dif®cult; at present, we cannot give a general recipe for the construction of scaling functions. Also, sampling must be ef®cient (an inherent problem of MD simulations) and this might limit the applicability of the method; nevertheless, the approach presented is promising and the assistance of more sophisticated sampling methods [31] might help to extend the investigations to more complex compounds. Acknowledgements We gratefully acknowledge support for B.N. from the Deutsche Forschungs-gemeinschaft under Grant no. DFG Kn 417/12 and general ®nancial support by EMS-Chemie AG, Domat/Ems, Switzerland. Generous amounts of computer time was provided by the Competence Centre for Computational Chemistry (C4) of the ETH ZuÈrich. We also acknowledge valuable discussions with Andrei A. Gusev. References [1] [2] [3] [4] [5] [6]
Shing KS, Gubbins KE. Mol Phys 1981;43:717. Takeuchi H, Okazaki K. J Chem Phys 1990;92:5643. Takeuchi H, Roe R-J, Mark JE. J Chem Phys 1990;93:9042. Mueller-Plathe F. Macromolecules 1991;24:6475. Gusev AA, Suter UW. Phys Rev A 1991;43:6488. Sok RM, Berendsen HJC, van Gunsteren WF. J Chem Phys 1992;96:4699. [7] de Pablo JJ, Laso M, Suter UW. J Chem Phys 1992;96:6157. [8] de Pablo JJ, Laso M, Suter UW. Macromolecules 1993;26:6180.
B. Nick, U.W. Suter / Computational and Theoretical Polymer Science 11 (2001) 49±55 [9] Pant PV, Boyd RH. J Chem Phys 1993;98:9895. [10] Gusev AA, Arizzi S, Suter UW, Moll DJ. J Chem Phys 1993;99: 2221. [11] Gusev AA, Suter UW. J Chem Phys 1993;99:2228. [12] Mueller-Plathe F, Rogers SC, van Gunsteren WF. J Chem Phys 1993;98:9895. [13] Han J, Boyd RH. Macromolecules 1994;27:5365. [14] Tamai Y, Tanaka H, Nakanishi K. Macromolecules 1994;27:4498. [15] Tamai Y, Tanaka H, Nakanishi K. Macromolecules 1995;28:2544. [16] Gee RH, Boyd RH. Polymer 1995;36:1435. [17] Gusev AA, Suter UW, Moll DJ. Macromolecules 1995;28:2582. [18] Ralston ARK, Denton DD, Bicerano J, Moll DJ. Comput Polym Sci 1996;6:15. [19] Bicerano J, Moll DJ. Comput Polym Sci 1996;6:117. [20] van der Vegt NFA, Briels WJ, Wessling M, Strathmann H. J Chem Phys 1997;105:8849.
55
[21] Knopp B, Suter UW, Gusev AA. Macromolecules 1997;30:6107± 6113. [22] Knopp B, Suter UW. Macromolecules 1997;30:6114±6119. [23] Insight 3.0.0. Discover 2.9.7. Amorphous Cell 8.0.0, MSI, San Diego, CA, 1996. [24] Verlet L. Phys Rev 1967;159:98. [25] Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. J Chem Phys 1984;81:3684. [26] Theodorou DN, Suter UW. Macromolecules 1985;18:1467. [27] Allen MP, Tildesley DJ. Computer simulation of liquids, Oxford University Press, 1987. [28] Brandrup J, Immergut EH, editors. Polymer handbook 3. New York: Wiley, 1989. [29] Flory PJ. J Chem Phys 1942;10:51. [30] Huggins ML. J Chem Phys 1941;9:440. [31] Leontidis E, Suter UW. Mol Phys 1994;83:489.