Free-energy landscape of polymer-crystal polymorphism

Polymorphism rationalizes how processing can control the final structure of a material. The rugged free-energy landscape and exceedingly slow kinetics in the solid state have so far hampered computational investigations. We report for the first time the free-energy landscape of a polymorphic crystalline polymer, syndiotactic polystyrene. Coarse-grained metadynamics simulations allow us to efficiently sample the landscape at large. The free-energy difference between the two main polymorphs, $\alpha$ and $\beta$, is further investigated by quantum-chemical calculations. The two methods are in line with experimental observations: they predict $\beta$ as the more stable polymorph at standard conditions. Critically, the free-energy landscape suggests how the $\alpha$ polymorph may lead to experimentally observed kinetic traps. The combination of multiscale modeling, enhanced sampling, and quantum-chemical calculations offers an appealing strategy to uncover complex free-energy landscapes with polymorphic behavior.


I. Introduction
The complex interplay of molecular interactions can lead to a material exhibiting multiple distinct forms in its solid state. 1 This polymorphism often results in wildly different material properties, making the study of polymorphism both essential for quality control in manufacture, and also a fascinating structure-property problem.Beyond structure and property, the intermediate processing of a material has a key impact on the resulting polymorph.Fundamentally, this stems from two ingredients: (i) the underlying free-energy landscape being sufficiently rugged to display several low-lying metastable states; and (ii) the exceedingly slow kinetics exhibited in the solid phase, preventing a full/ergodic kinetic relaxation.
The screening of polymorphs has traditionally exclusively been performed experimentally, in spite of the significant costs involved.Computational methods hold the promise of predicting polymorphic stability before going to the laboratory.][4][5][6][7][8][9][10][11][12][13] The modeling of polymorphism holds two challenges: sampling and modeling accuracy.The free-energy landscape exhibits an overwhelming number of configurations, of which only an infinitesimal fraction competes in terms of low-lying states.Furthermore, correctly ranking the relative free energies of each conformer requires computational methods that are accurate enough (of the order of thermal energy, k B T) to reproduce the underlying interactions.][16] Here, we focus on polymers, which not only embody countless industrial applications, but also for which we critically lack a detailed picture of their free-energy landscape.Evidently, the increased number of atoms per molecule involved will lead to higher structural correlations, significantly larger barriers, and much longer timescales compared to molecular crystals of small molecules.While many semicrystalline polymers possess only a single type of unit cell, there are exceptions: syndiotactic polystyrene (sPS), for instance, is well-known for its complex crystal polymorphism.
In this work, we aim at unraveling the underlying freeenergy landscape of sPS, to better understand the body of experimental evidence gathered around its polymorphs.We focus on the thermally-induced processing aspect -the a and b polymorphs, shown in Fig. 1. [18][19][20][21][22][23] They share the same intrachain conformations, but with different interchainpacking structures.The a and b forms of sPS are further classified into limiting disordered forms, a 0 and b 0 , and limiting ordered forms, a 00 and b 00 . 18,20,22,24The processing conditions impact the forms experimentally observed. 22,24Fig. 1 displays the limiting ordered forms.
5][26][27][28] They hint at a highly complex freeenergy landscape.The relative stability between the a and b polymorphs is noteworthy: (i) starting from a high initial temperature, fast (slow) cooling will lead to the a (b) forms; (ii) under an identical slow cooling rate, melt crystallization starting under 230 1C and above 260 1C will yield the a and b polymorphs, respectively, while intermediate temperatures generate mixtures thereof.These results suggest that given sufficient mobility thanks to a high initial temperature and a slow enough cooling rate, the preferred packing structure corresponds to the b form.In the case of stiffened chains and/or reduced molecular mobility, the a polymorph is preferred.These observations of the structural properties of sPS indicate that crystallization to the a form results from a kinetically-controlled process, while b would be the thermodynamically stable form.Our understanding, therefore, falls short in several ways: how the free-energy landscape translates into the apparent differences between the two forms, and also a more mechanistic insight regarding the origin of these effects.
][36] These studies helped understand the structural features of some of these forms.Unfortunately, the atomistic resolution involved strongly limits the timescale that can be reached with the simulations -on the order of nanoseconds.This prevents both the observation of self assembly and also polymorph interconversion, thereby hindering access to the free-energy landscape.
To address the time-scale issue, we turn to coarse-grained (CG) modeling.By lumping several atoms into one larger superparticle or bead, CG models can sample significantly faster, while offering a systematic connection to the reference chemistry. 37Some of us recently applied a structure-based CG model aimed at reproducing certain thermodynamic aspects of sPS. 17Despite a parametrization and validation performed exclusively in the melt, we found remarkable transferability to the crystalline phase: not only does the CG model stabilize the a and b polymorphs, the melting temperatures of the two phases were found to be in excellent agreement. 38Our study aimed at an exploration of the self-assembly mechanisms of sPS, using a temperature-based enhanced-sampling molecular dynamics (MD) technique -parallel tempering.In the present work, we instead turn to methods based on collective variables (CVs), specifically metadynamics. 39,40We will show that an appropriate choice of CVs can lead to convergence of the simulations, and offer us access to a great diversity of structural forms.We present for the first time the free-energy landscape of polymorphism of a polymer crystal.
We further challenge the calculations of the free-energy difference between a and b polymorph stability by means of quantum-chemical calculation at the density functional theory (DFT) level.The results show excellent agreement with the CG simulations given the change in resolution.Critically, we find in both cases the preferential stabilization of the b phase -in line with experiments.

A. Metadynamics
In the Methods section, we present collective variables (CVs) that are capable of distinguishing five different phases of sPS (Section IV B).A significant distinction between these phases is essential to also enable the discovery of other intermediate phases.A further requirement is the absence of hidden barriers that would hinder dynamics along the CVs. 39To alleviate possible artifacts due to an inappropriate choice of CVs, we tested several of them and later re-weighted all simulations to the same CV space.This further allows us to empirically check the convergence of our simulations.
We focus on a two-dimensional CV exploration, as a balance between exploration and convergence: a three-dimensional CV-space exploration can require excessive memory and presents challenges to converge due to the curse of dimensionality.][43][44] We herein present four combinations of CVs (referenced in Table S1, ESI †) based on orientational features between the v 2 transverse vectors (Fig. 1b).Variants of S focus on specific angles, while the second Legendre polynomial P 2 probes the overall ordering: Fig. 2 shows a metadynamics simulation at T = 400 K driven by the combination of CVs DS & P 2 ðv 2 Þ using two walkers.The walkers were initiated from the two crystalline phases a and b.The simulations are able to transition many times between the main polymorphs, starting at around 20 ns (Fig. 2a).The results highlight that the ability of the CVs to distinguish a from b helps ensure large conformational transitions (Fig. S2, ESI †).We note however the absence of conformations in the amorphous phase, due to both our sampling below the transition temperature (roughly 450 K) and our protocol's restraining box range and chain direction (see Section IV C).

B. Polymorphic stability
Having established that our combination of CVs enables a satisfactory transition between major polymorphs, we turn to the question of convergence.To assess convergence, we run metadynamics in the CV space of our four combinations, and marginalize the free-energy landscape to only display their common CV: the SMAC difference DS.Fig. 3 compares the free-energy surface as a function of the CV GðDSÞ.We have marked the a and b polymorphs according to values of the CV from unbiased MD simulations (see the ESI †).All four curves agree within roughly 5 kJ mol À1 across the range of DS values, despite their sampling along different complementary CVs.Convergence as a function of simulation time is further displayed in Fig. 4a, which focuses on the free-energy difference between the a and b polymorphs, G a À G b .We find that all curves converge after roughly 1 to 2 ms.We do see variations between simulations reminiscent of the spread in panel (a).
Given the remarkable complexity of probing the free-energy landscape of polymer crystals, we consider this level of agreement as an encouraging indicator of the level of convergence of our simulations.
Metadynamics simulations on different system sizes lead to similar free-energy profiles (see Fig. S5, ESI †), whether changing the number of chains in the box or the number of monomers per chain.We rationalize the lack of system-size   dependence in two ways: (i) the lack of dependence on the number of chains can be explained by the collective nature of the transitions, where the box is so small that all chains transition to a new phase at once; and (ii) the free-energy barriers do not scale with the number of monomers, because the transitions are orthogonal to the chain direction, as indicated by the order parameters P 2 (v 1 ) and P 2 (v 2 ) (see Fig. S2, ESI †).As such, the size regime we work with appears to have a lack of system-size dependency, despite the extensivity of the free energy.Identifying this scaling would require the simulation of much larger simulation boxes, for instance because multiple structures could occur (concurrent with grain boundaries).
Having identified the two polymorphs as local minima with a rationalization of kinetic routes in between them, we try to further establish their relative thermodynamic stability.For this, we start from the experimentally determined structures and employ quantum mechanical methods for local structure relaxations and prediction of their temperature dependent free energy (see Section IV D).In contrast to the molecular dynamics simulations, our DFT results do not describe anharmonicities of the energy surface, potentially neglecting some thermal effects.On the other hand, the described interactions are at a quantum-mechanical level, physically more sound, and thus expected to be more accurate compared to the classical potentials used in our molecular dynamics.We note that while the a and b polymorphs for the quantum-mechanical calculations are taken from experimental structures, the coarse-grained references stem from the simulations.Our recent analysis of the CG model indicated a faithful stabilization of the a structure but discrepancies in the b structure, likely due to sidechain packing issues. 38Using experimental and simulation b structures for the quantum-mechanical calculations and CG simulations, respectively, allows us to account for the limited resolution of the CG model.
We show in Fig. 4b the quantum-mechanical energy difference between the two polymorphic forms.In a static picture at 0 K, sPS b is predicted to be more stable by 4 kJ mol À1 .Upon heating to 400 K, the enthalpy difference increases by 1 kJ mol À1 , while the free energy difference decreases by 3 kJ mol À1 .At elevated temperatures of about 500 K, form a is predicted to be more stable, which is, however, not experimentally observable due to the amorphous phase.The stability difference at room temperature is in excellent agreement with the estimations from our dynamics simulations, indicating a stabilization of b by 5-10 kJ mol À1 .The smaller stability difference predicted by DFT is in line with our experience from molecular crystals, where higher quality interaction energy models typically lead to smaller energy gaps between polymorphs. 14The overall analysis matches the experimental expectation that form a crystallizes on rapid cooling from elevated temperatures, while b forms in slow cooling experiments. 45he analysis reported in Fig. 4 shows that the b form is systematically better stabilized than the competing a form.The 1D landscape clearly separates a from b at the left and right sides of the range, respectively.These are separated by both a/b mixtures and the amorphous phases at around DS % 0. Interestingly, we observe a significantly lower free-energy barrier upon going from the amorphous phase to the pure a polymorph, compared to that for the b polymorph: while the former is between 10 and 15 kJ mol À1 high, the other is up to 20 kJ mol À1 .The b form is more stable across the CV space, but the a form is easier to obtain from the mixture and amorphous phases -a kinetic effect.This can help rationalize the kinetic-trap behavior of the a form found experimentally. 45Some of us previously identified an overpopulation of the a form when probed in simulation box geometries concomitant to the a unit cell, suggesting a templating mechanism. 38

C. Free-energy landscape of sPS
As an extension to Fig. 3, Fig. 5 shows a representation of the free-energy landscape for the CV combination DS & S 1 .Stability is color coded from blue to red.We observe a large diversity of phases with distinct structural features.Notably, we find many more structures than our previous study based on parallel tempering. 38e first analyze structures similar to the a polymorph.Of particular interest are symmetries around triplets of chains that form the fundamental symmetric unit of a structures (see Fig. 1b).As compared to our previous work that used parallel tempering, we observe a broader variety of relative orientations between chain triplets.Small but noticeable variations can be found among the a-type structures on the landscape.The differentiation between a 0 and a 00 is made more difficult by imposing 12 chains in our simulation box, while the unit cell of a contains only 9 chains. 20We emphasize the difference between apparently distinct forms stabilized in our simulations: a 0 , a 0 defect and a 00 defect , shown in Fig. 5.All triplets of chains, represented by groups of tan-colored beads, exhibit virtually identical orientations in a 0 and little variation in a 0 defect : the angles between side-chain vectors are almost all strictly at 1201, leading to DS % À0:4. a 00 like , on the other hand, displays two different triplet orientations, analogous to the experimentally resolved a 00 .The angles between side-chain vectors do not always correspond to 1201, leading to DS % À0:2, located near the mixture phases.Our simulations do not stabilize the a 00 polymorph, which may be due to the number of chains being incongruent with the unit cell, or possibly limitations in the CG force field in reproducing fine steric features. 17,38n line with our previous study, we find an alternate form to the experimentally-resolved 22,24 limiting disordered b 0 and limiting ordered b 00 polymorphs as the most global minimum of sPS: b sim , where we highlighted the difference in layering. 38hile the parallel tempering simulations led to neither experimental form, the metadynamics simulations successfully sampled them, albeit at too high free energy: the b 0 exp and b 00 exp forms sit at about 30 and 50 kJ mol À1 higher than the global minimum, respectively.We argued before that the simple description of the side-chain sterics likely had a detrimental effect on the stability of the b polymorphs.This effect was motivated by a discrepancy in the melting temperature from CG simulations, as compared to reference atomistic simulations: in excellent agreement for the a form, but with too-low stability for b.In the context of the present work, these structural artifacts likely give rise to shifts in the free-energy landscape shown in Fig. 5.

III. Conclusion
In this paper, we study the free-energy landscape of syndiotactic polystyrene (sPS) using computational methods to get microscopic insight into polymorphic interconversion.The development of adequate collective variables (CVs) applied to metadynamics, together with the use of a remarkably transferable coarse-grained (CG) model, allows us for the first time to cross the significant barriers between polymorphs in polymer crystals.Minute structural differences between polymorphs require finely-tuned CVs to account for the small variations.Rather than relying on a single CV combination, running metadynamics simulations on several such combinations and re-weighting them help us test for convergence of the simulations.We find excellent agreement between four such combinations, despite the significant barriers exerted by the system.
We rely on a combination of two different SMAC variables 46 to build a free-energy landscape.The b form clearly stands as the global minimum, even though we observe fine differences between the simulated and experimental layerings, arguably an artifact of the side-chain representation in the CG model.Encouragingly, we do observe the two experimental b 0 exp and b 00 exp structures in the metadynamics runs.The a form is between 5 and 10 kJ mol À1 less stable than the global minimum.As a complementary approach, we used quantum mechanical models to compute the temperature dependent relative stability of the a and b polymorphs.We predict the stability to be slightly smaller (1 kJ mol À1 ), but can overall confirm the CG simulations.
Remarkably, we find a significantly lower free-energy barrier upon going from the amorphous phase to a rather than b.This lowered activation could partially explain the a polymorph's tendency to be identified as a kinetic trap.It also complements our previous observation: a is overstabilized in box geometries congruent with its unit cell, a templating mechanism of sorts.Differences in nucleation rates may further help drive the system in the direction of a, although this is beyond the scope of this work.
Varying the system size will naturally affect the results: smaller simulation boxes are more likely to suffer from incompatibilities with different crystal units, resulting in artificially low stabilization.On the other hand, the significant free-energy barriers will become more challenging to cross as simulation boxes grow.The present study demonstrates that a multiscale approach can provide insight and complementary information to experiments on polymer-crystal polymorphism.

A. Coarse-grained simulations
We rely on a previously-developed coarse-grained (CG) model for syndiotactic polystyrene (sPS), referred to hereafter as the Fritz model. 17It maps each monomer onto two types of CG beads: ''A'' for the chain backbone and ''B'' for the phenyl ring (Fig. 1a).The model represents PS as a linear chain of alternating A and B CG beads, supplemented by sophisticated bonded potentials to ensure accurate structures, including the correct tacticity -enabling the crystallization of sPS.The bonded interactions were obtained by direct Boltzmann inversion of distributions obtained from atomistic simulations of single chains in a vacuum.The nonbonded potentials are derived using the conditional reversible work (CRW) method: 47,48 constraineddynamics runs with the all-atom model of two short chains in a vacuum.

B. Collective variables
Metadynamics acts on the selected CVs to drive the system and cross free-energy barriers.This puts an essential emphasis on the choice of CVs, known to be extremely system and process dependent. 39][54] Crystallization is typically characterized by long-range order.For example, cluster symmetries are often probed via the Steinhardt parameter, Q l , 55,56 which represents the rotationallyinvariant spherical harmonics of order l.Q 4 and Q 6 , in particular, have been studied in the context of Lennard-Jones particles, 52 ice, 53 and calcium carbonate nanoparticles. 54To study polymorphism in sPS, however, we have found the Steinhardt parameter to inefficiently distinguish crystalline forms (see Fig. S3, ESI †).We rationalize this by the lack of differentiation for transverse vectors (see Fig. 1b).This has led us to the development of CVs that are tailored to sPS.
Fig. 1b shows a longitudinal view of a chain conformation, as well as transverse sections of the two main polymorphs of interest: the experimentally-determined a and b structures. 24e herein propose two sets of vectors: one longitudinal vector, v 1 , and one transverse vector, v 2 .v 1 is oriented along the backbone, and it is defined as the interparticle vector between two consecutive monomers (i.e., backbone beads).As for the transverse vector v 2 , for each monomer, it originates from the second backbone and points to the bisector of the two closest side chains.Fig. 1b indicates that typical angles for the a and b polymorphs are roughly 1201 and 1801, respectively.Note that here, to improve the efficiency of computation, the side-chain vector is an average over each whole molecule.Effectively this assumes a homogeneous configuration across the chain, in line with the system sizes considered here.
Because these two polymorphs consist of identical chain conformations, no differences can be extracted from intrachain statistics.Instead, differences between polymorphs are concentrated in the interchain configurations.As such, any CV that is able to distinguish between polymorphs ought to focus on interchain geometries.By symmetry, the most noticeable differences occur along the transverse sections.These differences can be displayed more intuitively by the distribution of angles between neighboring transverse vectors v 2 (see Fig. S1, ESI †), which is derived from unbiased MD simulations at 300 K.In the following discussion, we will introduce two kinds of CVs that are functions of these transverse angles.
1. Legendre polynomial P 2 .The second Legendre polynomial P 2 probes the orientation between two vectors e i and e j : P 2 is a natural candidate to describe orientational ordering and has been used to monitor the crystalline growth and/or state of polymer chains. 57,58Some of us previously used P 2 to monitor the melting transition of sPS. 38. SMAC.Giberti et al. 46 recently introduced a CV to capture the inherent variety of crystal symmetries. 59This CV, simply called SMAC in Plumed, was formulated with the aim of accounting for both local density and the mutual orientation of molecules.In our case, it probes features of the angle distributions between transverse vectors, as shown in Fig. S1 (ESI †).It compares an input angle, y ij , with a reference angle, y n , via K n (y ij À y n ) = e À((yijÀyn) 2 /2sn 2 ) . ( This kernel smoothly interpolates from identical to distant angles leading to values from 1 to 0, respectively.SMAC relies on these kernels to probe one or multiple reference angles according to the phase of interest, supplemented by a smooth cutoff scheme: Eqn (3) relies on the two switching functions sðrÞ ¼ and c(r) = exp(Àr/r c ), where r s and r c (see Table S1, ESI †) provide a balance to keep the CV local, while incorporating enough numbers.The quantity is then averaged over all s i : Based on the angle distributions of Fig. S1 (ESI †), we construct a number of SMAC CVs using various sets of reference angles to optimally distinguish between sPS phases.Table S1 (ESI †) lists the set of CVs we used in this work.For instance, the main peaks for the a and b phases led to the definition of SMAC CVs S a and S b , centered at 1201 and 1701, respectively.
To emphasize both features at once, we constructed a CV based on their difference: DS ¼ S b À S a .This can be observed by monitoring the CV during unbiased MD simulations (Fig. S2c, ESI †).Three additional SMAC CVs are presented in Table S1 (ESI †).

C. Metadynamics
Well-tempered metadynamics simulations were performed at 400 K. 40,[60][61][62][63][64][65] The bias deposition stride was set to be 0.5 ps, and the bias factor was 20.The Gaussian-bias width was set to 0.01 for DS, and 0.05 otherwise.The initial Gaussian heights were all set to 3 kJ mol À1 .Each simulation consisted of multiple walkers: 66 In metadynamics, we refer to a biased run that explores the CV-space as a walker, while multiple walkers simultaneously explore and carve the same free-energy landscape.This method can significantly speed up the convergence of the simulations, as all walkers contribute to a single, combined free-energy landscape.In this work, 2 walkers were run in parallel and initialized from the a and b forms, respectively.All simulations were performed using GROMACS 5.1.4 67and PLUMED 2.4. 68,69Simulations were carried out in the isothermal-isobaric ensemble at P = 1 bar using the velocity-rescale thermostat 70 and the anisotropic Parrinello-Rahman barostat. 71We ensured stable variations in the simulation box by restraining the range of allowed geometries (5 o a 2 o 9; 5 o b 2 o 9; 6 o c 2 o 8 nm 2 ).We also restrained the direction of each chain to lie within 301 of the box's z component.This avoided significant collective rotations of the chains with respect to the simulation box, leading to alignments along the other coordinates.Such a scheme merely aims at forcing all chains to loosely lie within an arbitrarily-chosen z axis.This helped avoid artifacts when calculating longitudinal and transverse vectors, especially at higher temperatures.More details can be found in the ESI.† Here, E latt (V) is the zero-temperature internal energy of the crystal given per monomer unit-the lattice energy.
where the phonon frequencies o k,p correspond to a k-point in the first Brillouin zone and a phonon band index p.The temperature-dependent harmonic enthalpy, H HA , is described in the ESI.3][74][75] The DFT calculations are done with a screened exchange hybrid density functional, dubbed HSE-3c. 76,77It combines accurate descriptions of geometries over a broad class of systems with an efficient treatment of non-local exchange and long-range London dispersion interaction. 78,79For a general overview of dispersion corrections in the density functional framework and the treatment of molecular crystals, see ref. 80-82.The implementation of HSE-3c into the CRYSTAL17 program enables the fast computation of electronic structures and phonon modes using all point-and space-group symmetries. 83,84Geometry optimizations are performed with tight convergence thresholds and in space groups P3 1 (a) and Pnma (b).The Brillouin zone has been sampled with a 1 Â 1 Â 5 and 1 Â 5 Â 3 grid for the a and b polymorphs, respectively.G-Point frequencies have been computed by sHF-3c, 85,86 and the above G-point frequencies are found to be negligible at the DFTB3-D3 level, [87][88][89] which has been shown to be a reliable approach for organic solids. 90

Fig. 1
Fig.1(a) Molecular structure of polystyrene and the CG-mapping scheme of the Fritz model.17(b) Left: Longitudinal view of an all-trans chain conformation; Right: transverse section of the experimentallyresolved a and b forms.In the transverse sections, each polymer chain is represented by three CG beads (highlighted in red dashed circles).The red solid lines represent cross-sectional vectors pointing from the backbone to the bisector of its two closest side chains.

Fig. 2
Fig. 2 Time evolution of the CVs during two metadynamics runs (each walker shown with a different color) using the combination of (b) DS & (c) P 2 (v 2 ).(a) First 60 ns of DS, represented by a gray area in panel b.

Fig. 3
Fig. 3 Convergence of the free-energy calculations.Comparison of four metadynamics simulations at T = 400 K with different CVs (see labels) projected on DS.

Fig. 4
Fig. 4 Free-energy difference between a and b forms.(a) The CG simulations show the time evolution of G a À G b from the different metadynamics simulations.(b) The DFT-based stability differences at standard pressure and as a function of temperature, including lattice energy E latt , enthalpy H HA , and free energy G HA .
D. Quantum-chemical methodsQuantum-chemical methods are used to simulate the relative thermodynamical stability of the two sPS polymorphs.The two polymorphic forms identified by the metadynamics simulations and experimentally characterized in the literature have been taken as starting points for local geometry optimization.Phonon modes are computed to confirm the stationary points as local minima and to give access to the temperaturedependent harmonic Gibbs free energy G HA (T,P) = E latt + G HA vib (T) + PV.