Introduction

Coarse-grained (CG) molecular dynamics (MD) simulations are increasingly playing an important role in our efforts to understand how macroscopic properties of glassy polymers depend on their molecular architecture. Such understanding takes us closer to the important technological goal of designing the architecture to achieve a set of targetted macroscopic properties.

Coarse graining drastically reduces the degrees of freedom compared to an all-atom MD simulation by lumping selected atoms into super-atoms or beads [1,2,3], which serve as interaction centres in computations. This has the potential to extend the time scale significantly and the length scale modestly when MD computations are conducted with the CG polymer. Coarse-graining techniques differ in the methods they use to calibrate the bonded and non-bonded force fields [2, 4,5,6,7,8,9]. All coarse-graining approaches are aimed at capturing specific static and dynamic properties of the underlying all-atom polymeric system [10,11,12,13,14,15]. CG MD simulations have enabled researchers to investigate problems such as polymer indentation [16], mechanics of crosslinked amorphous polymer adhesives [17], mechanical behaviour of polymer thin films [18], and thermomechanical properties of bulk polymers [19].

We have been interested in determining stress–strain responses of glassy, amorphous polymers starting from a description of their molecular architectures [20,21,22]. It has been shown that uniaxial stress–strain responses of glassy polymers (e.g. polystyrene and a class of polyimides) predicted from such carefully calibrated CG simulations do indeed qualitatively resemble their experimentally obtained responses. However, since uniaxial stress–strain responses of most glassy thermoset and thermoplastic polymers are sensitive to the applied strain rate and the CG simulations are conducted at strain rates that are many orders of magnitude higher than those used in experiments, the stress-carrying capacities of the polymers are greatly overestimated, even though the features of the stress–strain responses are reproduced faithfully.

From CG simulations mostly on polymer melts, it is known that dynamic properties obtained from CG simulations tend to differ from detailed all-atom simulations. This is because the smooth beads generally lead to faster dynamics in the CG system. As a result, the evolution of the detailed system and the all-atom one over time differ. Seen in another way, an unit of time in the CG system may correspond to a much longer time in the detailed system. Time- and temperature-scaling techniques [23,24,25,26,27] have been devised to utilize the faster dynamics to extend the time scales accessible by CG simulations. Similarly, techniques to extend the CG force fields calibrated at a particular temperature to other temperatures have also been proposed [28,29,30,31,32,33].

In this paper, we explore the following question:

While CG simulations conducted at extremely high rates cannot be expected to quantitatively predict the experimental stress–strain response of a glassy polymer, is it possible to predict qualitative differences between the responses of two polymers? For instance, given two molecular architectures and similar CG schemes, is it possible to predict which of the two will have higher strength?

To this end, we have chosen two polymers with similar molecular architectures as test cases. These are, glassy, amorphous polyether ether ketone (PEEK) and polyether ketone ketone (PEKK). Both of these are thermoplastic polymers which find applications in various fields, such as aerospace, manufacturing, automotive, etc. These two polymers have the same functional groups in their molecular structures, but in different proportions, as depicted in Fig. 1. The ratio of ether to ketone groups is 2:1 in PEEK, whereas it is 1:2 in PEKK. An early all-atom MD study of PEEK [34] suggested that the ether-ether rings are more mobile than the ether–ketone rings. Also, the ketone groups are more rigid than the ether groups. Thus, due to a greater proportion of ketone groups, PEKK chains are overall expected to be more rigid, and consequently PEKK has a glass transition temperature \({T}_{\mathrm{g}}\) that is about \(20\,^{\circ }\mathrm{C}\) higher than PEEK [35,36,37].

Figure 1
figure 1

Molecular structure of the monomer of (a) PEEK and (b) PEKK. The coarse-graining method adopted for the two polymers are also shown. The aromatic benzene and the ketone (C=O) groups have been considered as beads ‘A’ and ‘B’, respectively, whereas the ether group (–O–) has been considered as bead ‘O’.

Full size image

However, PEEK and PEKK have very similar mechanical properties on account of their similar molecular architectures. Maintaining mechanical properties while offering a superior thermal resistance, PEKK can potentially replace PEEK in many applications. Will a CG scheme that sufficiently reproduces static properties of PEEK and PEKK be able to capture the fact that both have similar uniaxial stress–strain response at a given strain rate?

Similar questions have been asked by other researchers in different contexts. For example, Tschöp et al. [2] coarse grained three varieties of polycarbonate melts and showed that (by extrapolating the high-temperature diffusion constant to the zero diffusion limit) the Vogel Fulcher temperatures are in line with experimental values. Two different CG mapping schemes for polystyrene involving different apportioning of atomic masses have been used by Harmandaris et al. [12] to demonstrate that the ratio of the masses of the beads has a bearing on the characteristic ratios and radial distribution functions (rdfs). As the number of beads used to coarse grain an all-atom polyimide is systematically increased, enhanced friction on individual beads leads to a better approximation of the uniaxial stress–strain response even without using artificially enhanced friction or time scaling [21]. The effect of the number of beads has been studied to determine both static and dynamic properties of polyethylene [27].

We have considered a coarse-graining scheme (see Fig. 1) where we assume three types of beads and six beads per monomer. This model considers the aromatic rings (type A beads), ether (type O), and ketone groups (type B) as different beads. In this case, the point of difference between PEEK and PEKK is the sequence of CG beads along the backbone. The reduction in the number of atoms per monomer from atomistic to CG model is from 34 to 6 in PEEK and from 35 to 6 in PEKK.

The remainder of this paper has been organized as follows. In Sect. “Results and discussion”, we critically compare various predictions from CG simulations on the two polymers. Salient conclusions are presented in Sect. “Conclusions”. The CG procedure adopted here is similar to many others and a very brief description is given in Sect. “Computational methods”. Special emphasis is given to the differences in force field calibrations between PEEK and PEKK. Owing to their similar architecture, the differences are subtle.

Results and discussion

Reported properties of PEEK and PEKK

In amorphous form, the reported densities of PEKK [38] and PEEK [39] are \(1270\) and \(1264 \mathrm{kg}/{\mathrm{m}}^{3}\) respectively. But the glass transition temperature of PEKK is \(160\,^{\circ } \mathrm{C}\) compared to \(147\,^{\circ } \mathrm{C}\) for PEEK.

The intrinsic uniaxial stress–strain response of amorphous glassy polymers can be obtained from compressive tests where various localization phenomena like necking and crazing do not intervene [40]. Compressive strength of PEEK is well characterized over a wide range of strain rates [41,42,43]. At room temperature, PEEK is rate sensitive in compression with the compressive strength ranging from about \(125 \mathrm{MPa}\) at \(\dot{\epsilon }=1{0}^{-4} {\mathrm{s}}^{-1}\) to \(175 \mathrm{MPa}\) at \(1{0}^{2} {\mathrm{s}}^{-1}\) [43]. It shows some hardening at large compressive strains only when the rate of loading is high.

On the other hand, we could not locate a systematic study of rate dependence of compressive behaviour of amorphous PEKK. Compressive strength of amorphous PEKK reported in commercial datasheets [38] is slightly lower than amorphous PEEK—about \(108 \mathrm{MPa}\) at a low strain rate. The rate dependence of the compressive response of PEKK has not been explored in the literature. The higher glass transition temperature, combined with its comparable compressive strength, makes PEKK an important substitute for PEEK in several applications.

Coarse-grained approximation of the structure of PEEK and PEKK

We start by examining the averaged end-to-end distance of the all-atom and CG molecules. The distance between the centre of mass of the terminal monomers is considered as the end-to-end distance in both all-atom and CG molecules. In view of the fact that PEKK contains a greater proportion of the stiff ketone (C=O) linkages, its end-to-end distance is expected to be more than PEEK. For the all-atom as well as the CG model, the end-to-end distance of PEKK does indeed turn out to be larger than PEEK. The end-to-end distance in the CG model is close to that of all-atom structures. The end-to-end distance in 100-monomer-long PEEK is \(\sim\) 58 \(\text{\AA }\) in all-atom, with a standard deviation of 10.9 \(\text{\AA }\), and \(\sim\) 59 \(\text{\AA }\) in CG chains, with a standard deviation of 10 \(\text{\AA }\). For 100-monomer-long PEKK, the end-to-end distance is \(\sim\) 69 \(\text{\AA }\) in all-atom, with a standard deviation of 9.5 \(\text{\AA }\), and \(\sim\) 66 \(\text{\AA }\) in CG chains, with a standard deviation of 14.5 \(\text{\AA }\) (Supplementary Information, Fig. 3). The six beads per monomer ratio used in the CG model seems sufficient to capture the long-range structure of the polymer chains.

The fidelity of the molecule obtained in the CG scheme also becomes evident when we plot the rdfs of the monomers of PEEK and PEKK, shown in Fig. 2. The rdf of the monomers in the CG sample match that of the all-atom case remarkably well for both PEEK and PEKK. The closeness of the end-to-end distances and monomeric rdfs computed from the CG model to those in the all-atom cases show that a 6 beads per monomer mapping of PEEK and PEKK is adequate, as far as approximation of the long-range structures of the molecules is concerned.

Figure 2
figure 2

Radial distribution functions (rdf) of monomers of (a) PEEK (solid lines) and (b) PEKK (dashed lines) in the CG model compared with those from the all-atom system (red curves).

Full size image

Properties that depend on rate

We have established that the CG systems are structurally very close to corresponding all-atom ones as well as to each other. Consequently, amorphous, glassy PEEK, and PEKK have structural properties that are almost identical. But, in spite of this similarity, predicting properties like glass transition temperature and uniaxial stress–strain response of amorphous, glassy systems through CG simulations require careful considerations.

Though the energies of the CG and all-atom systems are matched at the initial state, their subsequent dynamics is different and they may evolve differently over time. Consequently, evolution of rate-dependent physical quantities with time, in general, will not be the same in the all-atom and CG systems.

The variation of the density with temperature is a case in point. At room temperature, the densities of CG and all-atom systems of PEEK and PEKK are identical. This is seen in Fig. 3a and b. However, when the samples are heated to a higher temperature at a rate \(\dot{T}\), the densities start to diverge from each other. As shown in Fig. 3a and b, the change in density between room temperature and \(500\,^{\circ }\mathrm{C}\) is maximum for the all-atom systems. With coarse graining, the density becomes less sensitive to temperature changes.

Figure 3
figure 3

Variations of the density of a PEEK and b PEKK with temperature are shown for the all-atom and the CG model. The intersections of the straight lines fitted to the low-temperature and high-temperature data points give the value of \({T}_{\mathrm{g}}\), which are marked by vertical lines dropped onto the horizontal axes.

Full size image

Also, the glass transition temperatures in the all-atom systems of PEEK and PEKK can be obtained by the familiar construction shown in Fig. 3a and b. The intersection of the straight lines fitted to the low- and high-temperature data in the density versus temperature plots gives the value of \({T}_{\mathrm{g}}\). The value of \({T}_{\mathrm{g}}\) obtained for PEEK (Fig. 3a) is about \(30\,^{\circ } \mathrm{C}\) lower than that for PEKK (Fig. 3b). This is consistent with experiments [35, 44]. Expectedly, the absolute values of \({T}_{\mathrm{g}}\) obtained for both PEEK and PEKK are higher than the experimental values due to the very high rates of change of temperature used in the MD simulations [45].

As the decrease in density with temperature is gentler for CG systems, the change in slope at \({T}_{\mathrm{g}}\) is hardly visible. Similar results have been reported by other authors. Coarse graining makes the determination of \({T}_{\mathrm{g}}\) difficult. In case of polyimides, Pandiyan et al. [21] had to perform coarse graining with 8 beads per monomer in order to discern the value of \({T}_{\mathrm{g}}\). Significantly lower values of the coefficient of thermal expansion for CG polystyrene has also been reported by Hsu et al. [15]. In our case, the CG scheme does exhibit a barely discernable slope change, yielding values of \({T}_{\mathrm{g}}\) that are somewhat higher than for the all-atom system.

The example above demonstrates the difficulties in predicting rate-dependent properties with CG simulations. ‘Time’ in the CG system is different from the physical time in the all-atom one. Thus, the \(\dot{T}\) that the all-atom and CG samples ‘feel’ are different. Hence, the density versus temperature behaviour is also different in each case.

The effective rate of change of a physical quantity in a CG system can be estimated if the scaling between physical time and that in the CG system is known. The mean squared displacement (msd) is often used to determine the scaling between the physical and CG time. The msd is defined as

$$\langle {\Delta }^{2}(t)\rangle =\frac{1}{N}\sum_{i}|{\mathbf{r}}_{i}(t)-{\mathbf{r}}_{i}(0){|}^{2},$$
(1)

where \({\mathbf{r}}_{i}(t)\) denote the position vectors of the centre of mass of the \(i\)th chain, or the \(i\)th monomer in the sample. In polymer melts, the scaling between the times is often determined by demanding that the evolution of the msd of the CG chain centroids in scaled time is identical to that of the all-atom chain centroids in physical time.

Unfortunately, in glassy state, the msd of the chain centroids or, for that matter monomers, is extremely small and remains almost constant over time. This is apparent from Fig. 4 where the evolution of msd of the monomeric units in all-atom and CG systems are shown. It is clearly not possible to use the monomeric msd’s to determine the scaling between physical time and the time in the CG system. Alternate routes to determine the scaling need to be explored for glassy systems.

Figure 4
figure 4

The mean squared displacement of the monomers of PEEK (solid lines) and PEKK (dashed lines) in the all-atom and CG model is plotted in log–log scale against time.

Full size image

Let us turn to uniaxial compressive stress–strain responses of the all-atom and CG polymers. The responses of both PEEK and PEKK as well as the values of their compressive strengths are strongly strain rate-dependent. The response in uniaxial deformation can be described by

$$\sigma =C(\epsilon ,T){\left(\dot{\epsilon }\right)}^{m},$$
(2)

where \(\sigma\) and \(\epsilon\) are the uniaxial stress and strain, respectively. The strain rate exponent is denoted by \(m\). The prefactor \(C\) can be a function of strain \(\epsilon\) and temperature \(T\). It determines the shape of the uniaxial stress–strain curve, while the second term on the right decides the extent of hardening with strain rate.

Expectedly, the all-atom systems of PEEK and PEKK [shown in Fig. 5a and b respectively] exhibit almost similar uniaxial compressive responses. The compressive strength of both PEEK and PEKK is \(\sim 200 \mathrm{MPa}\), about twice the experimental value at low strain rate. The uniaxial compressive simulations are conducted in a \(N{\sigma }_{xx}{\sigma }_{yy}{L}_{z}\) ensemble, where the virial stresses \({\sigma }_{xx}={\sigma }_{yy}=0\) (simulating a pure uniaxial situation) and \({L}_{z}\) is reduced to half its initial value (i.e. a compressive strain of \(0.5\) is imparted) in \(0.1 \mathrm{ns}\). For the all-atom system, this translates to an extremely high strain rate of about \(5 {\mathrm{ns}}^{-1}\).

Figure 5
figure 5

Uniaxial compressive response of (a) PEEK and (b) PEKK at \(300 \mathrm{K}\) are shown. In both the figures, solid curves pertain to all-atom simulations. For PEEK in (a), the all-atom response is plotted at rate of \(5 {\mathrm{ns}}^{-1}\). For PEKK in (b) all-atom responses are shown at \(50\) and \(5 {\mathrm{ns}}^{-1}\). The dashed curves pertain to the CG scheme. For PEEK, these are shown for \(5\) and \(0.5 {\mathrm{ns}}^{-1}\). For PEKK, responses at \(5, 0.5, 0.05\) and \(0.005 {\mathrm{ns}}^{-1}\) are shown.

Full size image

The CG systems, when deformed at the same applied strain rate of \(\dot{\epsilon }=5 {\mathrm{ns}}^{-1}\), produce stress–strain responses that qualitatively resemble the all-atom as well as experimental ones. For example, the peak compressive stress occurs at a strain of about 10% as in the all-atom and a hardening response is seen at high compression.

However, quantitatively at \(\dot{\epsilon }\) = 5 \({\mathrm{ns}}^{-1}\), the compressive strengths for PEEK and PEKK in the CG systems are considerably higher than the corresponding all-atom systems. Also, though the all-atom systems of PEEK and PEKK exhibit similar compressive strengths (a fact seen in experiments, as discussed earlier), in the CG systems at the same applied strain rate of \(5 {\mathrm{ns}}^{-1}\), PEKK has higher strength than PEEK.

The scaling between the physical strain rate and the CG rate can be determined from relatively small-sized simulations on all-atom and CG samples. To this end, variations of the strengths (the maximum stress borne by a sample) with strain rate are plotted on a log–log scale in Fig. 6a–c. The plot for PEEK in Fig. 6a shows two slopes, indicating two different values of the exponent \(m\) in Eq. 2. The slopes denoting the values of \(m\) are indicated on the log–log plots and are approximately same for the all-atom and CG cases. In case of PEEK (see, Fig. 6a), to achieve a uniaxial response identical to that of an all-atom simulation at a physical strain rate of \(\dot{\epsilon }=5 {\mathrm{ns}}^{-1}\), the CG system will have to be deformed at a strain rate about 10 times slower. The scaling factor \(\alpha\) between the physical strain rate and that felt by the CG system is approximately \(10\).

Figure 6
figure 6

Variations of the maximum stress-carrying capacity or strength of (a) PEEK and (b) PEKK with uniaxial strain rate \(\dot{\epsilon }\) are plotted on a log–log scale. In (c) the variations for PEEK and PEKK are compared for the CG case. In (a) and (b) the slopes of the curves, i.e. the values of the strain hardening exponent \(m\), are indicated. The horizontal arrows show the approximate correspondence between a physical strain rate of \(\dot{\epsilon }=5 {\mathrm{ns}}^{-1}\) in the curve on the right to the effective strain rate on the left.

Full size image

Guided by the simple analysis above, the CG system for PEEK is run at a strain rate one order of magnitude lower than the all-atom. As shown in Fig. 5a, this brings the entire uniaxial compressive stress–strain response of PEEK very close to that of the corresponding all-atom system simulated at \(5 {\mathrm{ns}}^{-1}\).

The value of \(m\) is lower than PEKK, indicating a somewhat weaker strain rate hardening. As seen from Fig. 6b, to mimic the response of an all-atom sample at \(\dot{\epsilon }=5 {\mathrm{ns}}^{-1}\), the CG sample has to be subjected to an effective strain rate much smaller than \(0.005\), the slowest rate at which we could conduct the simulations. But, an all-atom sample deformed at \(\dot{\epsilon }=50 {\mathrm{ns}}^{-1}\) corresponds to a CG sample deformed at \(0.005 {\mathrm{ns}}^{-1}\). This is again shown in Fig. 5b, where the solid curve for \(\dot{\epsilon }=50 {\mathrm{ns}}^{-1}\) matches the dashed one for \(0.005 {\mathrm{ns}}^{-1}\).

The fact that PEEK and PEKK have almost similar uniaxial responses cannot be seen from CG simulations if both are subjected to the same strain rate. As seen from Fig. 6c, a CG sample of PEEK subjected to \(5 {\mathrm{ns}}^{-1}\) should correspond to a CG sample of PEKK subjected to a rate between \(0.5\) and \(0.05 {\mathrm{ns}}^{-1}\). As is apparent from Fig. 5a and b, the entire uniaxial response in Fig. 5a at \(5 {\mathrm{ns}}^{-1}\) falls between curves pertaining to \(0.5\) and \(0.05 {\mathrm{ns}}^{-1}\) in Fig. 5b.

Due to closeness of the molecular architectures of PEEK and PEKK, a closer correspondence between strain rates was expected between their CG simulations in Fig. 6c. However, as indicated by the arrow in Fig. 6c, the factor is close to 1000. The strain rate sensitivity of PEKK, to the best of our knowledge, has not been characterized experimentally. However, the large difference in strain rate sensitivities of PEEK and PEKK seems to be an artefact of the coarse-graining procedure. We believe that the strain rate sensitivity of such CG systems is related to the rate at which LJ non-bonded interactions break and reform when the ensemble is strained.

A final point about the CG scheme is in place. The force fields for this scheme were calibrated at \(27\,^{\circ } \mathrm{C}\), while the glass transition temperature is determined as \(\sim 250\,(280)\,^{\circ } \mathrm{C}\) in PEEK (PEKK). In the glassy regime, the forcefields seem to be transferable as far as predicting the uniaxial response is concerned. This is shown in the Supplementary Information.

Conclusions

CG MD simulations of long chained amorphous, glassy polymers can be used for predicting their mechanical properties. The question that we ask in this work is whether, given two polymers with similar architectures, a suitable CG model can be designed that allows us to make a comparative assessment of their mechanical properties? This question is important because often, like in the case of PEEK and PEKK that we studied, subtle changes in the structure of the monomer are seen to lead to technologically significant changes in key properties. It is useful to be able to predict if a proposed change in the molecular architecture will lead to, e.g. higher glass transition temperature or mechanical strength.

The first step to answer the above question is designing a CG scheme that reproduces the overall structure of the all-atom system with high fidelity. We have shown that for PEEK and PEKK, markers of the long-range structure like monomeric rdf’s and average end-to-end lengths are reproduced well with a 6 beads per monomer coarse-graining scheme with 3 different beads.

Though the structure and initial energy of the glassy polymers are approximated well, rate dependence of many physical properties like glass transition temperature and mechanical strength poses further challenges. This is mainly a result of the extremely sluggish dynamics of dense, glassy systems. In such cases, quantitative prediction of strongly rate-dependent properties like uniaxial stress–strain response is uncertain. These challenges arise from the fact that time in CG systems do not represent physical time. The scaling between the physical time and that in CG systems depend on the degree of coarse graining as well as the molecular architecture of the polymers. Even slightly different architectures like PEEK and PEKK may lead to significantly different scalings (see, Figs. 5 and 6).

As far as uniaxial stress–strain response is concerned, in order to obtain the scaling between a physical strain rate in an all-atom system and apparent strain rate in the corresponding CG one, we propose conducting a set of all-atom and CG simulations at different strain rates (see, Eq. 2). From these simulations, the factor \(\alpha\) that decides the effective strain rate \(\dot{\epsilon }/\alpha\) that a CG sample should be subjected to in order to mimic uniaxial response at a physical strain rate of \(\dot{\epsilon }\) can be determined (see, Fig. 6).

Our proposed coarse-graining scheme, while providing computational advantages, leads to good approximation of the average end-to-end distance of the polymer chains and the general structure. The rdfs are closely approximated and even dynamic properties with weak rate sensitivity like glass transition temperature can be predicted.

It is, however, clear that CG systems at glassy temperatures find it difficult to capture the rate sensitivity of their uniaxial stress–strain behaviour. As future work, one of the following two modifications are worth persuing:

  1. (1)

    An iterative Boltzmann inversion technique for the non-bonded interactions, leading to a tabulated non-bonded potential (rather than the harder LJ potential used by us). However, it must be noted that, the time in CG systems is still expected to be different from real time. A time scaling between all-atom and CG will still be necessary.

  2. (2)

    Use of time-dependent LJ parameters as used by Hsu et al. [15] and Song et al. [31] could help in capturing the rate dependence of uniaxial response better.

Computational methods

The procedures adopted for calibrating CG force fields for PEEK and PEKK are similar to those reported earlier [2, 4, 11, 22, 46]. All MD simulations have been conducted with the large-scale atomic/molecular massively parallel simulator (LAMMPS) [47]. In cases where temperature or pressure in the ensembles are constant, Nose–Hoover thermostat or barostat has been utilized. Velocity Verlet algorithm has been employed to integrate the Newton’s equations of motion. The steps involved in the process are discussed in Sects. “Step 1: All-atom simulations with PEEK and PEKK” through “Comparison between the calibrated CG forcefields of PEEK and PEKK” with a time step, \(\Delta\)t, of 0.1 and 5 fs for all-atom and coarse-grained simulations, respectively. The CG model simulations are 3–4 times faster than the all-atom simulations due to fewer interaction centres and combined with the 10–50 times increase in timestep, the CG simulations offer a speedup of around 50–200, i.e. O(100). The following subsections describe the preparation method of the all-atom and coarse-grained samples.

Step 1: All-atom simulations with PEEK and PEKK

Physical quantities obtained from detailed all-atom MD simulations on PEEK and PEKK serve as benchmarks that the CG simulations reported later should match up to. The procedures adopted towards this end, for the two polymers, are exactly identical.

  1. (1)

    Initially, an energy-minimized chain of PEEK and PEKK, consisting of 100 monomers is prepared using the “Amorphous Cell” module of Biovia Materials Studio, employing the polymer consistent force field, PCFF [48].

  2. (2)

    Inside a periodic box, 10 copies of this chain are placed. The box contains 10 chains, each 100 monomers long with a total of 3402 atoms in PEEK and 3502 in PEKK.

  3. (3)

    MD simulations using the NPT ensemble are performed on LAMMPS [47]. After initial equilibration at \(300 \mathrm{K}\) and \(1 \mathrm{atm}\) pressure, the starting density of the sample turns out to be quite low. To bring its density closer to the experimentally observed density, the sample is taken to a very high pressure (\(\sim 1{0}^{4} \mathrm{atm}\)) in steps to increase the density to a value sufficiently higher than the experimentally reported density of each polymer.

  4. (4)

    The pressure is then reduced to one atmosphere. If the desired density is not attained, the pressure is again increased to a higher value. The cycling between the high pressure and atmospheric is continued till the experimentally observed density under atmospheric conditions is attained. In all-atom samples, attaining close to the correct experimentally observed density at glassy temperatures and atmospheric pressure is often difficult [49]. A similar strategy has been used by Lyulin et al. [50] to prepare all-atom samples of polyimides.

  5. (5)

    Once the desired density is achieved, the sample is further equilibrated in an NPT ensemble at \(300\mathrm{ K}\) and one atmosphere pressure. The all-atom equilibration simulations have been performed for \(\sim 1000\mathrm{ ps}\). The energy evolution plots for PEEK and PEKK are provided in Fig. 1 of Supplementary Information. During this process, the density and total energy remained unchanged while the bond lengths, angles and dihedrals fluctuated slightly about their targetted equilibrium values. During the entire sample preparation procedure, the chains moved distances in the order of their lengths, causing change in density of two orders of magnitude (from \(0.04\) to the desired \(1.2 \mathrm{g}/{\mathrm{cm}}^{3}\)).

Step 2: Coarse graining: calibrating the bonded force field

The CG configuration is characterized by a set of bonded parameters \(R,\varTheta\) and \(\varPhi\). In the CG scheme, the bond length \(R\) can denote \({R}_{\mathrm{OA}}\) or \({R}_{\mathrm{AB}}\). The bond angles \(\varTheta\) may be \({\varTheta }_{\mathrm{OAO}},{\varTheta }_{\mathrm{BAB}},{\varTheta }_{\mathrm{OAB}},\dots\) while the dihedrals \(\varPhi\) represents \({\varPhi }_{\mathrm{AOAO}},{\varPhi }_{\mathrm{OABA}},{\varPhi }_{\mathrm{AOAB}}\dots\).

The basic assumption inherent in calibrating the bonded potentials is that the scaled probability distributions for the bond length \({P}_{\mathrm{S}}(R, T)\), angle \({P}_{\varTheta }(\varTheta , T)\) and dihedral \({P}_{\varPhi }(\varPhi , T)\) at a temperature \(T\) are independent so that the joint probability distribution factorizes as [2]:

$$P(R,\varTheta ,\varPhi ,T)={P}_{\mathrm{S}}(R, T){P}_{\varTheta }(\varTheta , T){P}_{\varPhi }(\varPhi , T).$$
(3)

The validity of this assumption is often verified once the CG procedure is complete [11, 12, 21]. The scaled probability distributions are obtained from unscaled ones as [11],

$${P}_{\mathrm{S}}(R, T)\propto \frac{{{P}_{\mathrm{S}}}^{\prime}(R, T)}{{R}^{2}},$$
(4)

for bond lengths,

$${P}_{\varTheta }(\varTheta , T)\propto \frac{{{P}_{\varTheta }}^{\prime}(\varTheta , T)}{\mathrm{sin}\varTheta }$$
(5)

for bond angles, and,

$${P}_{\varPhi }(\varPhi , T)\propto {{P}_{\varPhi }}^{\prime}(\varPhi , T),$$
(6)

for dihedrals. The unscaled probability distributions are denoted by a prime. For a generic configuration variable \(\zeta\) (where \(\zeta =R,\varTheta\) or \(\varPhi\)), the unscaled probability distributions are determined by simply counting \(N(\zeta )\), the number of instances where the variable \(\zeta\) lies between \(\zeta -\Delta\) and \(\zeta +\Delta\), \(\Delta\) being sufficiently small.

A simple Boltzmann inversion [2, 51] then yields the respective free energy differences, i.e.

$$V(\zeta )\propto -{k}_{\mathrm{B}}T\mathrm{ln}p(\zeta ).$$
(7)

Here, \({k}_{\mathrm{B}}\) is the Boltzmann constant and \(V(\zeta )\) is a measure of the free energy difference, or potential energy in vacuum [4]. Finally, the bonded potentials are fitted with class 2 forcefield equations provided in LAMMPS [47].

Step 3: Coarse graining: calibrating the non-bonded force field

To calibrate the non-bonded force fields, beads are placed at the geometric centres of the group of atoms in the all-atom sample which the bead represents. This implies that 6000 beads of 3 types (A, B or O) are placed. The beads A, B and O have masses 76, 28 and 16 respectively.

Consider a bead pair formed by species \(I\) and \(J\) (where \(I\) and \(J\) can be one of A, O and B). We denote a generic \(IJ\) pair by \(\alpha\). Consider that the ensemble has \({N}_{\alpha }\) such bead pairs,Footnote 1 with the distance between the centroids of the \({{i}_{\alpha }}\mathrm{th}\) pair denoted by \({R}_{{i}_{\alpha }}({i}_{\alpha }\in [1,{N}_{\alpha }])\). Then the non-bonded interaction between the beads \(I\) and \(J\) is taken to have the form

$$U\left({\sigma }_{\alpha }, {\epsilon }_{\alpha };{ R}_{{i}_{\alpha }}\right)=\left\{\begin{array}{ll}{\epsilon }_{\alpha }\left[2{\left(\frac{{\sigma }_{\alpha }}{{R}_{{i}_{\alpha }}}\right)}^{9}-3{\left(\frac{{\sigma }_{\alpha }}{{R}_{{i}_{\alpha }}}\right)}^{6}\right]& \mathrm{for}\,{R}_{{i}_{\alpha }}< {R}_{\mathrm{c}} \\ 0& {R}_{{i}_{\alpha }} \ge {R}_{\mathrm{c}},\end{array}\right.$$
(8)

where \({R}_{\mathrm{c}}\) is the cut-off distance.

If \(p\) atoms of the all-atom polymer are mapped onto bead of type \(I\) and \(q\) in \(J\), we fit the parameters \({\sigma }_{\alpha }\) and \({\epsilon }_{\alpha }\) by minimizing, in a least squares sense,

$$\Pi \left({\sigma }_{\alpha }, {\epsilon }_{\alpha }\right)=\sum_{{i}_{\alpha }=1}^{{N}_{\alpha }}{\left[U\left({\sigma }_{\alpha }, {\epsilon }_{\alpha }; {R}_{{i}_{\alpha }}\right)-\sum_{i=1,i\in I}^{p}\sum_{j=1,j\in J}^{q}{U}_{ij}^{\mathrm{PCFF}}\left({r}_{ij},{\sigma }_{ij},{\epsilon }_{ij}\right)\right]}^{2}.$$
(9)

Here \({U}_{ij}^{\mathrm{PCFF}}\) is the non-bonded potential in PCFF. So, the parameters \({\sigma }_{ij}, { \epsilon }_{ij}\) are known.

In practice, the above ‘energy matching technique’ is applied to cases when \(I=J\). Figure 7(a) and (b) depict the non-bonded energy matching of the benzene (type ‘A’) and ketone (type ‘B’) beads, respectively. When \(I\ne J\), and the pair \(IJ\) is designated by \(\alpha\), \({\sigma }_{\alpha }\) and \({\epsilon }_{\alpha }\) are determined from:

Figure 7
figure 7

Non-bonded interaction energy based on Eqs. 8 and 9 for (a) benzene bead (\(\alpha\) represents AA) and (b) CO bead (\(\alpha\) represents BB) in PEEK.

Full size image
$${\sigma }_{\alpha }={\left(\frac{{{\sigma }_{\beta }}^{6}+{{\sigma }_{\gamma }}^{6}}{2}\right)}^{1/6}$$
(10)

and,

$${\epsilon }_{\alpha }=\frac{2\sqrt{{\epsilon }_{\beta }{\epsilon }_{\gamma }}{{\sigma }_{\beta }}^{3}{{\sigma }_{\gamma }}^{3}}{{{\sigma }_{\beta }}^{6}+{{\sigma }_{\gamma }}^{6}}.$$
(11)

In the above, the pair \(II\) is denoted by \(\beta\) and the pair \(JJ\) is denoted by \(\gamma\). In Fig. 7a and b, the blue points are obtained from \({N}_{\alpha }\) pairs of type \(\mathrm{A}-\mathrm{A}\) or type \(\mathrm{B}-\mathrm{B}\) (i.e. interaction between two benzene rings or two ketone groups), each plotted at the distance (i.e. the distance between the centres of the benzene rings or ketone groups) at which it was detected in the equilibrated all-atom ensemble. The red curves are the fits assuming LJ interaction between the corresponding \(\mathrm{A}-\mathrm{A}\) or \(\mathrm{B}-\mathrm{B}\) beads.

We have heuristically modified the parameters in the CG LJ potentials to ensure that the target density is attained in the CG systems at an isotropic stress–state corresponding to atmospheric pressure. A discussion regarding the necessity of heuristic modification after energy matching is provided in the Supplementary Information. The parameters that affect density are \({\sigma }_{\mathrm{AA}}\) and \({\sigma }_{\mathrm{BB}}\), for benzene beads (type ‘A’) and ketone beads (type ‘B’) respectively. Contours of density of PEEK and PEKK on the \({\sigma }_{\mathrm{AA}}\) and \({\sigma }_{\mathrm{BB}}\) plane are shown in Fig. 8. The paths from the initial values obtained from the energy matching technique and the final values at which the desired density is attained, are shown with red lines.

Figure 8
figure 8

Contours of density of the CG system of (a) PEEK and (b) PEKK are plotted on the \({\sigma }_{\mathrm{AA}}-{\sigma }_{\mathrm{BB}}\) plane. The paths that we have chosen to reach the desired density from the one obtained after energy matching, are shown in red lines.

Full size image

All bonded and non-bonded calibrated parameters for CG model are presented in Table 1. The values of \(\sigma\) obtained for benzene beads are close to those obtained by Harmandaris et al. [12] (\(\sigma\)=5.2 \(\text{\AA }\)) and Hsu et al. [15] (\(\sigma\)=5.3–5.7 \(\text{\AA }\) for a temperature range of 200–500 K) for polystyrene molecules. Please note that atomic charges have not been considered in the CG model and it is assumed that the effect of atomic charges of the atomistic model is embedded in the non-bonded interaction between beads. Also, the monomer of PEEK and PEKK is charge neutral.

TABLE 1 Calibrated parameters for bonded and non-bonded forcefields for PEEK and PEKK pertaining to the CG model with three beads A, B, and O.
Full size table

Comparison between the calibrated CG forcefields of PEEK and PEKK

The bonded and non-bonded potentials obtained in the CG model, lead to stiff molecules. The bond AB, shown in Fig. 9a has a fixed length of \(3.3 \text{\AA }\) for both the molecules. Though not shown here, potential functions for bond AO is also identical for PEEK and PEKK. The angle potentials are also stiff. The potentials for the angle ABA that occurs in both the polymers is shown in Fig. 9b. Even though they are derived from separate coarse-graining procedures (as a result the parameters pertaining to the angle are slightly different for PEEK and PEKK), for both polymers the equilibrium value is \(\sim 10{0}^{\circ }\) for this angle. The potentials for angle OAO, which is present in PEEK and BAB, which appears only in PEKK are also shown in Fig. 9b. While the latter is slightly more flexible, both these angles have values of \(16{6}^{\circ }\). The non-bonded potential for beads of type ‘A’ is identical for both polymers (Fig. 9d). For beads of type ‘B’, the non-bonded potential has the same value of \(\sigma\) but PEKK has a somewhat deeper potential well.

Figure 9
figure 9

Potential functions calibrated at \(T=300 \mathrm{K}\) and normalized with \({k}_{\mathrm{B}}T\) are shown for the CG model. Solid lines pertain to PEEK and dashed to PEKK. The potentials shown are for (a) stretching of bond AB in both PEEK and PEKK, (b) bending of angles ABA in both polymers, OAO in PEEK and BAB in PEKK, (c) torsion of AOAO in PEEK and ABAB in PEKK and (d) non-bonded interactions for beads of type ‘A’ and ‘B’ in both

Full size image

The torsional potentials AOAB and OABA, that occur in both molecules, exhibit potentials that are identical in spite of the difference in parameters shown in Table 1. The difference between PEEK and PEKK manifests only in the torsion potentials for AOAO in PEEK and ABAB in PEKK (Fig. 9c). For PEEK, AOAO is a freely rotating dihedral angle. For PEKK, the angle ABAB is considerably stiffer with an equilibrium value close to \({0}^{\circ }\).

In summary, the difference in potential between PEEK and PEKK is limited to the dihedrals that are present in one but not in the other and a deeper potential well for beads of type ‘B’ in PEKK.