Developing a Transferable Coarse-Grained Model for the Prediction of

Mar 28, 2019 - In the present work, we develop a coarse-grained (CG) model for polyimide (PI) at 800K and 1atm by applying iterative Boltzmann inversi...
0 downloads 0 Views 1MB Size
Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Computational Chemistry

Developing a Transferable Coarse-Grained Model for the Prediction of Thermodynamic, Structural and Mechanical Properties of Polyimides at Different Thermodynamic State Points Chenchen Hu, Teng Lu, and Hongxia Guo J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00887 • Publication Date (Web): 28 Mar 2019 Downloaded from http://pubs.acs.org on March 28, 2019

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

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

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

Journal of Chemical Information and Modeling

Developing a Transferable Coarse-Grained Model for the Prediction of Thermodynamic, Structural and Mechanical Properties of Polyimides at Different Thermodynamic State Points Chenchen Hua,b, Teng Lua and Hongxia Guoa,b* a

Beijing National Laboratory for Molecular Sciences, Joint Laboratory of Polymer Sciences and

Materials, State Key Laboratory of Polymer Physics and Chemistry, Institute of Chemistry, Chinese Academy of Sciences, Beijing 100190, China b University

of Chinese Academy of Sciences, Beijing 100049, China

Email: [email protected] Abstract: In the present work, we develop a coarse-grained (CG) model for polyimide (PI) at 800K and 1atm by applying iterative Boltzmann inversion (IBI) and density correction method to derive the bonded and non-bonded interaction potentials. Although the CG force field is built at a single thermodynamic state point without any temperature correction, the CG model possesses a rather favorable temperature transferability in a wide temperature range of 300-800K at P=1atm and good pressure transferability to some extent in a certain pressure range from 0.1MPa to 30MPa. In addition to the local conformation and local packing distribution functions, the thermodynamic properties such as the glass transition temperature and the coefficient of linear thermal expansion are predicted correctly by the CG model, and the isothermal compressibility coefficients calculated from both atomic and CG models are on the same order of magnitude. Additionally, the stress-strain behavior under compression or tension of the CG model shows a qualitative agreement with the atomistic results and the corresponding values of elastic modulus of CG model at different temperatures roughly match with those of atomistic model. 1. Introduction Polyimides (PIs) are a kind of very important organic aromatic heterocyclic polymer with an imide group. They have been widely used as electrical insulators, substrates for flexible printed circuits, as well as insulation layers for semiconductor, due to their excellent heat resistance, mechanical, dielectric properties and so on.1-2 Polyimides are 1

ACS Paragon Plus Environment

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

prepared by polymerization of dianhydride and diamine monomers, and hence the properties of PIs can be controlled by tuning the chemical structures of dianhydride and diamine monomers and introducing some specific chemical groups3-6. Although the research of PIs has received considerable attention, the relationships between macroscopic properties and microstructures of PIs have not yet been understood clearly, which affects the molecular design and synthesis of novel polyimides and further hinders the improvement of material performance. For example, the role of chemical details or the stacking arrangement of molecular chains in the material properties of PI films such as the thermal expansion coefficient, the glass transition temperature, the mechanical properties and so on still remains unclarified.7 With the development of computer power and simulation algorithms, computer simulation techniques have evolved as a powerful tool to complement experimental evidence, predict the trends of macroscopic properties and establish the structure–property relations for polymeric systems8-14. To get insight into how the chemical structure of dianhydride and diamine monomers influences the performances of resulting polyimides, atomic-scale molecular dynamics (MD) simulations have been utilized to explore the relevant properties of PIs15-23. For example, Chi et al. 23 have adopted all-atom MD simulation to investigate the influence of local changes in polyimide structure on the performance of gas separation PI membranes (such as gas permeability, solubility and diffusivity) at the molecular level. However, when we focus on the mechanical behavior and viscoelastic properties of the materials, non-equilibrium MD simulations need to be carried out in relatively large systems on a long-time scale to obtain accurate statistical results, which is very computationally demanding by using atomistic MD simulations. In order to improve computational efficiency and reduce computation time, some studies of atomic-scale MD simulations often choose to turn off the electrostatic interaction of the system, which more or less affects the ability of model to quantitatively predict the material properties of PIs.22, 24 Therefore, it is crucial to establish an efficient computational model of polyimides that takes account of all interaction between atoms and allows one to access large spatial-temporal scales at a reasonable computational cost and hence 2

ACS Paragon Plus Environment

Page 2 of 43

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

Journal of Chemical Information and Modeling

quantitatively evaluate the physical and mechanical properties of system . Recently, the rapidly developing systematic Coarse-Grained (CG) methods provide feasible ways to construct the above model. In the creation process of a CG model, a group of chemically bonded atoms are coarse-grained into one CG bead and the interaction potentials or CG force fields (FF) are derived by applying “top-down” or “bottom-up” strategy. The systematic CG models not only effectively reduce the system’s degree of freedom, but also reasonably maintain the interaction between the chemical groups in atomistic MD simulations. Thus, the CG simulation is a promising tool to accurately predict the polymer structures and properties over a larger range of length and time scales with excellent computational efficiency. At present, CG method has been successfully employed to study the thermodynamic and structural properties of polymers such as polybutadiene25-28, polystyrene29-33, polyethylene34 and so on. So far, many useful coarse-graining methods have been devised for obtaining the CG force fields35-39 and can be classified into the following categories based on the target properties of research interest: structure-based method40-42 which aims to reproduce the structures generated by the atomistic MD simulations, force-based means36, 43-46 which matches the pairwise force distributions of CG beads with the original atomistic force distributions, and thermodynamic quantities based approaches13, 47-52 which reproduce the thermodynamic quantities (density, enthalpy, etc.) measured from experiments or atomistic MD simulations. Although the CG models constructed by the methods mentioned above can reproduce the target properties very well, their capacity for accurately predicting the non-target properties is fairly limited. Previous research has shown that CG potentials derived by traditional structure-based methods, which only include the two-body interaction but not the three- and multi-body interactions, can reproduce the structure obtained in atomistic MD simulation but cannot precisely predict the thermodynamic and mechanical properties of the system42, 53. For instance, Rahimi et al.

53

constructed a CG model of polystyrene through the structure-based

method and found that the Young's modulus measured by the CG model is nearly 20 times lower than that obtained from the atomistic simulation. The similar problem exists with the traditional force-based coarse-grained method. Rosch et al.54 developed the 3

ACS Paragon Plus Environment

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

CG model of polystyrene by employing the force-based CG approach and found that the stress-strain curve obtained by the CG model is far different from the atomistic results even though the force in the CG system matches well with that in atomistic system. With a further correction on the non-bonded potential, although the stress-strain curve from the CG model system is close to that from the atomistic simulation, the structure cannot be reproduced well. As we know, constructing a CG model with a wide spectrum of properties matching well with an underlying atomistic model is essential for the establishment of the structure-property relationship for polymers.25-26 This is especially true for the polyimide system, in addition to accurately reproducing the structural properties, the resulting CG model is required to predict the thermodynamic (the coefficient of thermal expansion (CTE) and the glass transition temperature (Tg)) and mechanical properties of PI as precisely as possible, due to their importance to evaluate the feasibility of PIs in engineering applications, such as adhesives for aircraft structures and turbo fan engine components55. At present, improving the representability (i.e., the ability to reproduce other properties which are not used in the parameterization of the CG force field) of CG models is one of the important challenges in the modeling community. Some studies have confirmed that combined techniques which include several target properties for the CG FF parameter fitting can more or less promote the representability of the CG model25,

56-58.

For example, by employing a novel force-based and structure-based

combined CG method, the resulting CG models of polystyrene and hexahydro-1,3,5trinitro-s-triazine can reproduce the mechanical properties (i.e., elastic modulus) as well as the molecular local conformations as the underlying atomistic results54. Additionally, Rossi et al.

33

have proposed thermodynamic-structural hybrid CG approach by

optimizing the MARTINI13 with the radius of gyration. The derived CG model of polystyrene can reproduce the thermodynamic (i.e., density) and relevant structural (i.e., structure factor) properties simultaneously. We note that a CG model which can replicate structures and thermodynamics obtained from atomistic simulations is very important to bridge levels of resolution in a multiscale simulation, wherein one has to switch between atomistic and CG levels of resolution or pass information between them 4

ACS Paragon Plus Environment

Page 4 of 43

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

Journal of Chemical Information and Modeling

in order to determine the overall properties of polymer materials.40 In addition, the recent work of our group further demonstrates that the combined structure-based and thermodynamic quantities-based CG method indeed has significant advantages in improving the representability of the CG model 25-26, 29. For instance, the systematic CG model of trans-1,4-polybutadiene25-26 which has been developed by the CG method as mentioned above, can reproduce not only the target properties including density and local structure (RDF and local conformation distribution functions), but also the nontarget properties such as the chain size, isothermal compressibility coefficient, bulk modulus, the storage modulus and loss modulus at the parameterized temperature. In light of these facts, it is quite valuable to construct a CG model of PI system by employing the thermodynamic-structural combined CG approach and examine whether it can predict structures, the thermodynamic quantities (especially Tg and CTE), and the mechanical

properties

(particularly

bulk

modulus

and

Young’s

modulus)

simultaneously. Research of this kind will not only help us to understand the link between structures and material properties but also provide insights into the molecular design of PIs for particular applications. As we know, to reproduce the structure of the underlying atomistic model system, the CG force field is commonly derived by using “bottom-up” strategy, which requires to carry out atomistic MD simulation at one thermodynamic state point. On the other hand, to predict some thermodynamic properties such as Tg and CTE as the atomistic model, the CG force field is required to possess good temperature transferability, as the values of Tg and CTE obtained in simulation are generally computed from the densitytemperature curve. Additionally, in order to reproduce the bulk modulus or another desirable thermodynamic property the isothermal compressibility coefficient by the CG model as accurately as possible, the CG force field should have pressure transferability to ensure that the response of bulk density to pressure is nearly the same as that of atomistic model. In fact, deriving CG potentials which exhibit both structural and thermodynamic consistency with the underlying atomic model over a wide range of thermodynamic condition will significantly advance the predictive capability of CG models. In principle, the CG force field usually exhibits a state-point dependency. 5

ACS Paragon Plus Environment

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

Page 6 of 43

However, when the difference between the weights of the removed degrees of freedom at different state points is not obvious, the CG potentials still possess state-point transferability.

54

For example, Moore et al.59 have constructed the CG model of n-

dodecane via multistate iterative Boltzmann inversion, which is suited to simulate systems over a range of thermodynamic states. Additionally, recent studies suggest the combined CG methods have advantages in improving the transferability of CG model.33, 52, 59-61

An-Tsung et al.62 combine the Sinoda-DeVane-Klein and IBI approaches to

derive the CG force field of perfluorosulfonic acid polymer membranes and the resulting CG model shows good temperature transferability. Similarly, Rossi et al.33 construct the CG model of polystyrene by thermodynamic-structure hybrid CG approach and the CG force field possesses a good transferability in a wide range of temperature 350~600K. Furthermore, the current work of our group also demonstrates that a combined structure-based and thermodynamic quantity-based CG method is beneficial to the improvement of the transferability of CG model and the CG models of trans-1,4-polybutadiene

and

polystyrene

can

reproduce

the

structural

and

thermodynamic properties obtained from the atomistic simulations at different thermodynamic state points.25-26,

29

For example, the CG model of trans-1,4-

polybutadiene reproduces the isothermal compressibility coefficient or bulk modulus within a certain level of accuracy and the CG model of polystyrene replicates Tg and CTE, indicating that the CG model constructed by combined CG method possesses transferability in a larger pressure and temperature range. Considering the advantages of thermodynamics–structure combined CG approach, it is worth to derive the CG force field of PI by such method at a single thermodynamic state point, wherein a good thermodynamic state (e.g. temperature and pressure) transferability is expected to be achieved to some extent. As addresses above, due to their unique functionalities, there has been considerable interest in predicting the material properties of polyimides based on their chemical structure. Among various polyimides, the advanced material PMDA/ODA, which is synthesized by pyromellitic dianhydride (PMDA) and 4,4′-oxydianiline (ODA), is a representative polyimide with good thermal and chemical stability and has been used 6

ACS Paragon Plus Environment

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

Journal of Chemical Information and Modeling

commercially since the early 1960s.63 The experimental data of PMDA/ODA are commonly available and conveniently taken as benchmarks for the estimation of atomistic model and CG force field. Therefore, we devise a systematic CG model of PMDA/ODA by using the structural and thermodynamic combined CG approach and examine whether it possesses good transferability and representability and can reproduce the thermodynamic, structural and mechanical properties obtained from atomistic MD simulations and the experiments at different thermodynamic state points or not. It is noteworthy that the establishment of reliable CG model of PMDA/ODA system can provide important guidelines and valuable references for the construction of other polyimide CG models. For example, we can develop the CG models of other novel polyimides by introducing some specific chemical groups or nanoparticles based on the CG model of PMDA/ODA and then predict the relevant physical and mechanical properties with computational effort orders of magnitude lesser than atomistic simulation, which can then inform the design and synthesis of new PIs or the design and fabrication of advanced PI composites according to the practical application. However, we should note that unlike the currently well studied polymers such as polystyrene and polybutadiene, in which the CG models are usually built by one or two types of coarse-grained beads due to the relatively simple chemical structures, we have to choose four types of CG beads to describe the PMDA/ODA in order to keep certain chemical details and guarantee computational efficiency. Thus, there are 10 nonbonded interactions between CG beads, which makes the process of CG FF parameterization rather complicated if using Lennard–Jones-type potentials to describe the non-bonded potentials although they can successfully describe the non-bonded CG bead–bead interactions in the CG models of polystyrene29, 64 and polybutadiene25-26. Therefore, we employ the iterative Boltzmann inversion (IBI) approach combined with density-correction (in a NPT ensemble) to obtain the non-bonded CG potentials in a tabulated form, and take both the radial distribution function (RDF) and the bulk density from the atomistic simulation at a given thermodynamic state point of T=800K and P=1atm as target properties during the CG process. The choice of such thermodynamic condition ensures the PI samples for the CG FF parameterization being in the melt state 7

ACS Paragon Plus Environment

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

and fast equilibrated, which speedups the optimization of CG force field. Besides, the similar method like the IBI approach combined with pressure-correction (in a NVT ensemble) has been applied to derive the CG non-bonded potentials of polymers such as polystyrene melts32 and oligomeric polyurea65, which can reproduce the bulk density at a certain temperature range. Accordingly, the present work is not only an experiment to exercise a CG method and estimate the capacity and limitations of the resulting CG model, but also it can provide insights into how to develop transferable and representative CG models for complex polymer systems so that they can be used in multiscale simulation studies of structure-property relations. The remainder of this paper is organized as follows. We describe the simulation method, CG model, and important quantities to characterize the properties of PMDA/ODA in section 2. Then in section 3, we present the simulation results and compare the results from CG simulations to those from atomistic MD simulations and experiments. Finally, we summarize the prime conclusions in the last section.

2. Methodology and simulation details 2.1 All-atom MD Simulation All-atom (AA) MD simulations are performed to generate the corresponding thermodynamic properties and structural distributions of PMDA/ODA systems as targets for the construction of CG model. The basic repeat unit of the PMDA/ODA

8

ACS Paragon Plus Environment

Page 8 of 43

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

Journal of Chemical Information and Modeling

(a)

(b)

N

C1

C2

N

O

C2

(c)

Fig. 1 (a) Chemical structure of the basic repeat unit of PMDA/ODA polyimide and coarse-grained mapping, and (b) figurative representation of a repeat unit of PMDA/ODA composed of coarse-grained beads. (c) The snapshot of atomistic model system of PMDA/ODA at T=800K and P=1atm. The oxygen, carbon, nitrogen, hydrogen atoms are shown in red, grey, blue and white, respectively. polyimide is illustrated in Fig 1(a). In our preliminary work, we cool down five atomistic model systems with different chain lengths of 4, 12, 20, 24, 28 monomer units by MD simulations with the same cooling rate and determine the glass transition temperatures of these five atomistic model systems from the inflection point of the hightemperature and low-temperature density data in the density-temperature curves. The resulting Tg values as a function of the monomer units per chain for the atomistic model system is shown in Fig S14. It can be seen that Tg for the system consisting of chains with a length of 12 monomer units is very close to in the platform region, i.e., the region where Tg is independent of molecular weight. This indicates that the chain length of 12 9

ACS Paragon Plus Environment

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

monomer units is sufficient that the end effect is negligible. In other words, we can use the sample with chain length of 12 monomer units and confidently make comparisons with experimental systems. Therefore, considering computational efficiency as well as the above fact, we packed 27 PMDA/ODA chains of 12 monomers per chain into an amorphous cell as the system to derive the CG force field, and the total number of atoms per AA sample is 12690. In our AA simulation, the box is sufficiently large with the length nearly two times bigger than average radius of gyration for PMDA/ODA chains. For example, at T=800K the value of radius of gyration is 19.370.06Å, whereas the box size is 54.950.09Å. Additionally, as typically demonstrated by the snapshot of atomistic model system Fig1 (c), the simulation box is large enough that a given molecule interacts with at most one image of a neighbor molecule. The effect of the periodic boundary conditions will not play a significant role. In fact, the choice of the atomistic system size is usually dictated by the two factors: on the one hand, the chain number should be sufficient to avoid the situation wherein one chain can interact at the same time with more than one periodic image of another chain; on the other hand, it should not be so many to incur too high computational efforts. The COMPASS force field66 which has been successfully used to simulate the properties of polyimides67-71 is employed in this work. The cut-off distances for the non-bonded interactions and the electrostatic interactions are both set to 2.0 nm, and the particle-mesh Ewald method is adopted to handle the long-range electrostatic interactions. Most AA MD simulations are performed in NPT ensemble by the Forcite module in Materials Studio (MS) software, wherein the Berendsen thermostat and barostat are used to keep the systems at constant temperature and pressure with coupling times of 0.1 and 0.5ps, respectively, and the time step is fixed to 1fs. We note that Ahunbay et al. performed NPT MD simulations with COMPASS force field and Berendsen thermostat/Barostat and studied the correlations between the properties (d-spacing, glass transition temperature, fractional free volume, etc.) and plasticization of fluorinated polyimide induced by CO2 sorption.17 Considering the similarity of simulated systems in our and Ahunbay’s work, we employ the same force field and Berendsen thermostat/Barostat as his. Note again that the atomistic model at T = 800K and P = 1 atm is used to derive the CG force field. 10

ACS Paragon Plus Environment

Page 10 of 43

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

Journal of Chemical Information and Modeling

Hence the parameterized temperature is larger than Tg, which is conducive to the rapid establishment of the CG force field. In order to evaluate the temperature transferability of CG force field, we need to compute the density-temperature curve to obtain the thermodynamic data (Tg and CTE) and characterize the distributions of structure to analyze the structural properties of PMDA/ODA from the atomistic simulations as benchmarks. In addition, the data can also be used to test the reliability of the AA force field. We choose 1.0×1011K/s as cooling rate. The range of cooling rate used for the atomistic system consists of 27 polyimide R-BAPS chains of 8 monomer units each is 1.5 × 1011~1.5 × 1013K/min, namely 2.5 × 109~2.5 × 1011K/s, within which our employed cooling rate for the atomistic model system of γ= 1.0×1011K/s is located. Lyulin et al.22 studied the relationship between the cooling rate and the glass transition temperature of polyimide R-BAPS and the resulting Tg at an infinitely slow cooling rate is in good agreement between with the experimental and computational value. They also found that the slope of the density-temperature curves in the low-temperature domain (T