Population Balance Equations for Reactive Separation in Polymer Upcycling.

Many polymer upcycling efforts aim to convert plastic waste into high-value liquid hydrocarbons. However, the subsequent cleavage of middle distillates to light gases can be problematic. The reactor often contains a vapor phase (light gases and middle distillates) and a liquid phase (molten polymers and waxes with a suspended or dissolved catalyst). Because the catalyst resides in the liquid phase, middle distillates that partition into the vapor phase are protected against further cleavage into light gases. In this paper, we consider a simple reactive separation strategy, in which a gas outflow removes the volatile products as they form. We combine vapor-liquid equilibrium models and population balance equations (PBEs) to describe polymer upcycling in a two-phase semibatch reactor. The results suggest that the temperature, headspace volume, and flow rate of the reactor can be used to tune selectivity toward the middle distillates, in addition to the molecular mechanism of catalysis. We anticipate that two-phase reactor models will be important in many polymer upcycling processes and that reactive separation strategies will provide ways to boost the yield of the desired products in these cases.

A population balance model for the kinetics of covalent organic framework synthesis.

This study presents a population balance model for the kinetics of nucleation and growth in covalent organic framework (COF) synthesis. The model incorporates second-order nucleation and first-order growth rates, consistent with proposals in the literature. Despite having non-linear terms, an implicit analytic solution is derived and then converted to explicit solutions for the monomer concentration and size distribution of COF flakes as a function of time. For experimental definitions of the induction time and the initial growth rate based on yield (y) vs time (t) curves, the model predicts power-law relationships: tind=0.409kN-1/3kG-2/3cA0-1 and dy/dtmax=0.965kN1/3kG2/3cA0, respectively. We discuss the implications for the interpretation of Arrhenius plots. We also discuss key discrepancies with experiments, including the predicted attainment of 100% yield instead of 30%-40% as observed and the value of the yield at the inflection point in the yield vs time curve. We suggest extensions to the model, including nucleation and growth kinetics with equilibrium solubility limitations and two-dimensional nucleation for the formation of multilayer COF particles.

Zirconium-Catalyzed C-H Alumination of Polyolefins, Paraffins, and Methane.

C-H/Et-Al exchange in zirconium-catalyzed reactions of saturated hydrocarbons and AlEt3 affords versatile organoaluminum compounds and ethane. The grafting of commercially available Zr(OtBu)4 on silica/alumina gives monopodal ≡SiO-Zr(OtBu)3 surface pre-catalyst sites that are activated in situ by ligand exchange with AlEt3. The catalytic C-H alumination of dodecane at 150 °C followed by quenching in air affords n-dodecanol as the major product, revealing selectivity for methyl group activation. Shorter hydrocarbon or alcohol products were not detected under these conditions. Catalytic reactions of cyclooctane and AlEt3, however, afford ring-opened products, indicating that C-C bond cleavage occurs readily in methyl group-free reactants. This selectivity for methyl group alumination enables the C-H alumination of polyethylenes, polypropylene, polystyrene, and poly-α-olefin oils without significant chain deconstruction. In addition, the smallest hydrocarbon, methane, undergoes selective mono-alumination under solvent-free catalytic conditions, providing a direct route to Al-Me species.

Role of stacking disorder in ice nucleation.

The freezing of water affects the processes that determine Earth's climate. Therefore, accurate weather and climate forecasts hinge on good predictions of ice nucleation rates. Such rate predictions are based on extrapolations using classical nucleation theory, which assumes that the structure of nanometre-sized ice crystallites corresponds to that of hexagonal ice, the thermodynamically stable form of bulk ice. However, simulations with various water models find that ice nucleated and grown under atmospheric temperatures is at all sizes stacking-disordered, consisting of random sequences of cubic and hexagonal ice layers. This implies that stacking-disordered ice crystallites either are more stable than hexagonal ice crystallites or form because of non-equilibrium dynamical effects. Both scenarios challenge central tenets of classical nucleation theory. Here we use rare-event sampling and free energy calculations with the mW water model to show that the entropy of mixing cubic and hexagonal layers makes stacking-disordered ice the stable phase for crystallites up to a size of at least 100,000 molecules. We find that stacking-disordered critical crystallites at 230 kelvin are about 14 kilojoules per mole of crystallite more stable than hexagonal crystallites, making their ice nucleation rates more than three orders of magnitude higher than predicted by classical nucleation theory. This effect on nucleation rates is temperature dependent, being the most pronounced at the warmest conditions, and should affect the modelling of cloud formation and ice particle numbers, which are very sensitive to the temperature dependence of ice nucleation rates. We conclude that classical nucleation theory needs to be corrected to include the dependence of the crystallization driving force on the size of the ice crystallite when interpreting and extrapolating ice nucleation rates from experimental laboratory conditions to the temperatures that occur in clouds.

Broken bond models, magic-sized clusters, and nucleation theory in nanoparticle synthesis.

Magic clusters are metastable faceted nanoparticles that are thought to be important and, sometimes, observable intermediates in the nucleation of certain faceted crystallites. This work develops a broken bond model for spheres with a face-centered-cubic packing that form tetrahedral magic clusters. With just one bond strength parameter, statistical thermodynamics yield a chemical potential driving force, an interfacial free energy, and free energy vs magic cluster size. These properties exactly correspond to those from a previous model by Mule et al. [J. Am. Chem. Soc. 143, 2037 (2021)]. Interestingly, a Tolman length emerges (for both models) when the interfacial area, density, and volume are treated consistently. To describe the kinetic barriers between magic cluster sizes, Mule et al. invoked an energy parameter to penalize the two-dimensional nucleation and growth of new layers in each facet of the tetrahedra. According to the broken bond model, barriers between magic clusters are insignificant without the additional edge energy penalty. We estimate the overall nucleation rate without predicting the rates of formation for intermediate magic clusters by using the Becker-Döring equations. Our results provide a blueprint for constructing free energy models and rate theories for nucleation via magic clusters starting from only atomic-scale interactions and geometric considerations.

Free energy barriers for anti-freeze protein engulfment in ice: Effects of supercooling, footprint size, and spatial separation.

Anti-freeze proteins (AFPs) protect organisms at freezing conditions by attaching to the ice surface and arresting its growth. Each adsorbed AFP locally pins the ice surface, resulting in a metastable dimple for which the interfacial forces counteract the driving force for growth. As supercooling increases, these metastable dimples become deeper, until metastability is lost in an engulfment event where the ice irreversibly swallows the AFP. Engulfment resembles nucleation in some respects, and this paper develops a model for the "critical profile" and free energy barrier for the engulfment process. Specifically, we variationally optimize the ice-water interface and estimate the free energy barrier as a function of the supercooling, the AFP footprint size, and the distance to neighboring AFPs on the ice surface. Finally, we use symbolic regression to derive a simple closed-form expression for the free energy barrier as a function of two physically interpretable, dimensionless parameters.

Size-Controlled Nanoparticles Embedded in a Mesoporous Architecture Leading to Efficient and Selective Hydrogenolysis of Polyolefins.

A catalytic architecture, comprising a mesoporous silica shell surrounding platinum nanoparticles (NPs) supported on a solid silica sphere (mSiO2/Pt-X/SiO2; X is the mean NP diameter), catalyzes hydrogenolysis of melt-phase polyethylene (PE) into a narrow C23-centered distribution of hydrocarbons in high yield using very low Pt loadings (∼10-5 g Pt/g PE). During catalysis, a polymer chain enters a pore and contacts a Pt NP where the C-C bond cleavage occurs and then the smaller fragment exits the pore. mSiO2/Pt/SiO2 resists sintering or leaching of Pt and provides high yields of liquids; however, many structural and chemical effects on catalysis are not yet resolved. Here, we report the effects of Pt NP size on activity and selectivity in PE hydrogenolysis. Time-dependent conversion and yields and a lumped kinetics model based on the competitive adsorption of long vs short chains reveal that the activity of catalytic material is highest with the smallest NPs, consistent with a structure-sensitive reaction. Remarkably, the three mSiO2/Pt-X/SiO2 catalysts give equivalent selectivity. We propose that mesoscale pores in the catalytic architecture template the C23-centered distribution, whereas the active Pt sites influence the carbon-carbon bond cleavage rate. This conclusion provides a framework for catalyst design by separating the C-C bond cleavage activity at catalytic sites from selectivity for chain lengths of the products influenced by the structure of the catalytic architecture. The increased activity, selectivity, efficiency, and lifetime obtained using this architecture highlight the benefits of localized and confined environments for isolated catalytic particles under condensed-phase reaction conditions.

Diffusion Attachment Model for Long Helical Antifreeze Proteins to Ice.

Some of the most potent antifreeze proteins (AFPs) are approximately rigid helical structures that bind with one side in contact with the ice surface at specific orientations. These AFPs take random orientations in solution; however, most orientations become sterically inaccessible as the AFP approaches the ice surface. The effect of these inaccessible orientations on the rate of adsorption of AFP to ice has never been explored. Here, we present a diffusion-controlled theory of adsorption kinetics that accounts for these orientational restrictions to predict a rate constant for adsorption (kon, in m/s) as a function of the length and width of the AFP molecules. We find that kon decreases with length and diameter of the AFP and is almost proportional to the inverse of the area of the binding surface. We demonstrate that the restricted orientations create an entropic barrier to AFP adsorption, which we compute to be approximately 7 kBT for most AFPs and up to 9 kBT for Maxi, the largest known AFP. We compare the entropic resistance 1/kon to resistances for diffusion through boundary layers and across typical distances in the extracellular matrix and find that these entropic and diffusion resistances could become comparable in the small confined spaces of biological environments.

Kinetic coefficient for ice-water interface from simulated non-equilibrium relaxation at coexistence.

In the theory of solidification, the kinetic coefficient multiplies the local supercooling to give the solid-liquid interface velocity. The same coefficient should drive interface migration at the coexistence temperature in proportion to a curvature force. This work computes the ice-water kinetic coefficient from molecular simulations starting from a sinusoidal ice-water interface at the coexistence temperature. We apply this method to the basal and prismatic ice planes and compare results to previous estimates from equilibrium correlation functions and simulations at controlled supercooling.

Oriented attachment kinetics for rod-like particles at a flat surface: Buffon's needle at the nanoscale.

The adsorption of large rod-like molecules or crystallites on a flat crystal face, similar to Buffon's needle, requires the rods to "land," with their binding sites in precise orientational alignment with matching sites on the surface. An example is provided by long, helical antifreeze proteins (AFPs), which bind at specific facets and orientations on the ice surface. The alignment constraint for adsorption, in combination with the loss in orientational freedom as the molecule diffuses toward the surface, results in an entropic barrier that hinders the adsorption. Prior kinetic models do not factor in the complete geometry of the molecule, nor explicitly enforce orientational constraints for adsorption. Here, we develop a diffusion-controlled adsorption theory for AFP molecules binding at specific orientations to flat ice surfaces. We formulate the diffusion equation with relevant boundary conditions and present analytical solutions to the attachment rate constant. The resulting rate constant is a function of the length and aspect ratio of the AFP, the distance threshold associated with binding, and solvent conditions such as temperature and viscosity. These results and methods of calculation may also be useful for predicting the kinetics of crystal growth through oriented attachment.

Free energies of crystals computed using Einstein crystal with fixed center of mass and differing spring constants.

Free energies of crystals computed using a center of mass constraint require a finite-size correction, as shown in previous work by Polson et al. [J. Chem. Phys. 112, 5339-5342 (2000)]. Their reference system is an Einstein crystal with equal spring constants. In this paper, we extend the work of Polson et al. [J. Chem. Phys. 112, 5339-5342 (2000)] to the case of differing spring constants. The generalization is convenient for constraining the center of mass in crystals with atoms of differing masses, and it helps to optimize the free energy calculations. To test the theory, we compare the free energies of LiI and NaCl crystals from calculations with differing spring constants to those computed using equal spring constants. Using these center of mass finite size corrections, we compute the true free energies of these crystals for different system sizes to eliminate the intrinsic finite-size effects. These calculations help demonstrate the size of these finite-size corrections relative to other contributions to the absolute free energy of the crystals.

Ion dissolution mechanism and kinetics at kink sites on NaCl surfaces.

Desolvation barriers are present for solute-solvent exchange events, such as ligand binding to an enzyme active site, during protein folding, and at battery electrodes. For solution-grown crystals, desolvation at kink sites can be the rate-limiting step for growth. However, desolvation and the associated kinetic barriers are poorly understood. In this work, we use rare-event simulation techniques to investigate attachment/detachment events at kink sites of a NaCl crystal in water. We elucidate the desolvation mechanism and present an optimized reaction coordinate, which involves one solute collective variable and one solvent collective variable. The attachment/detachment pathways for Na+ and Cl- are qualitatively similar, with quantitative differences that we attribute to different ion sizes and solvent coordination. The attachment barriers primarily result from kink site desolvation, while detachment barriers largely result from breaking ion-crystal bonds. We compute ion detachment rates from kink sites and compare with results from an independent study. We anticipate that the reaction coordinate and desolvation mechanism identified in this work may be applicable to other alkali halides.

Performing solvation free energy calculations in LAMMPS using the decoupling approach.

The decoupling approach to solvation free energy calculations requires scaling the interactions between the solute and the solution with all intramolecular interactions preserved. This paper reports a new procedure that makes it possible to these calculations in LAMMPS. The procedure is tested against built-in GROMACS capabilities. The model compounds chosen to test our methodology are ethanol and biphenyl. The LAMMPS and GROMACS results obtained are in good agreement with each other. This work should help perform solvation free energy calculations in LAMMPS and/or other molecular dynamics software having no built-in functions to implement the decoupling approach.

Tandem Catalysts for Polyethylene Upcycling: A Simple Kinetic Model.

Of all plastics, the most abundantly produced is polyethylene, most of which is destined for landfills, shipping ports, and natural environments. The limited degradability and recyclability of this synthetic polymer motivates the development of chemical recycling methods. One possible approach consists of selective depolymerization to propylene with tandem olefin metathesis and double bond isomerization catalysts. In this paper, we transform thousands of coupled rate equations, pseudo-steady-state approximations, and local density approximations into one simple and analytically solvable Fokker-Planck type equation. The Fokker-Planck equation gives concise expressions for the rate of propylene production and polymer molecular weight evolution as functions of catalyst concentrations, rate constants, and ethylene concentrations.

Absolute chemical potentials for complex molecules in fluid phases: A centroid reference for predicting phase equilibria.

Solid-fluid phase equilibria are difficult to predict in simulations because bound degrees of freedom in the crystal phase must be converted to free translations and rotations in the fluid phase. Here, we avoid the solid-to-fluid transformation step by starting with chemical potentials for two reference systems, one for the fluid phase and one for the solid phase. For the solid, we start from the Einstein crystal and transform to the fully interacting molecular crystal. For the fluid phase, we introduce a new reference system, the "centroid," and then transform to gas phase molecules. We illustrate the new calculations by predicting the sublimation vapor pressure of succinic acid in the temperature range of 300 K-350 K.

Solid-solid phase equilibria in the NaCl-KCl system.

Solid solutions, structurally ordered but compositionally disordered mixtures, can form for salts, metals, and even organic compounds. The NaCl-KCl system forms a solid solution at all compositions between 657 °C and 505 °C. Below a critical temperature of 505 °C, the system exhibits a miscibility gap with coexisting Na-rich and K-rich rocksalt phases. We calculate the phase diagram in this region using the semi-grand canonical Widom method, which averages over virtual particle transmutations. We verify our results by comparison with free energies calculated from thermodynamic integration and extrapolate the location of the critical point. Our calculations reproduce the experimental phase diagram remarkably well and illustrate how solid-solid equilibria and chemical potentials, including those at metastable conditions, can be computed for materials that form solid solutions.

How fluxional reactants limit the accuracy/efficiency of infrequent metadynamics.

In an infrequent metadynamics (iMetaD) simulation, a well-tempered metadynamics bias accumulates in the reactant basin, accelerating escapes to the product state. Like the earlier hyperdynamics strategy, iMetaD enables estimates of the unbiased escape rates. However, iMetaD applies the bias to visited locations in a collective variable (CV) space, not to the more specific visited locations in a full configuration space as done in hyperdynamics. This difference makes rate estimates from iMetaD sensitive to the choice of CVs, to parameters that control the bias deposition rate, and to the preparation of the initial state within the reactant basin. This paper uses an extremely simple discrete state model to illustrate complications that can arise in systems that exhibit fluxional transitions between sub-basins of the reactant state. Specifically, we show how the reactant-to-product escape time and relaxation times within the reactant basin(s) impose bounds on the admissible parameter choices for an iMetaD calculation. Predictions from the discrete state model are validated by iMetaD simulations on a corresponding two-dimensional potential energy surface.

Importance learning estimator for the site-averaged turnover frequency of a disordered solid catalyst.

For disordered catalysts such as atomically dispersed "single-atom" metals on amorphous silica, the active sites inherit different properties from their quenched-disordered local environments. The observed kinetics are site-averages, typically dominated by a small fraction of highly active sites. Standard sampling methods require expensive ab initio calculations at an intractable number of sites to converge on the site-averaged kinetics. We present a new method that efficiently estimates the site-averaged turnover frequency (TOF). The new estimator uses the same importance learning algorithm [Vandervelden et al., React. Chem. Eng. 5, 77 (2020)] that we previously used to compute the site-averaged activation energy. We demonstrate the method by computing the site-averaged TOF for a simple disordered lattice model of an amorphous catalyst. The results show that with the importance learning algorithm, the site-averaged TOF and activation energy can now be obtained concurrently with orders of magnitude reduction in required ab initio calculations.

Solvent reaction coordinate for an SN2 reaction.

We study the prototypical SN2 reaction Cl- + CH3Cl → CH3Cl + Cl- in water using quantum mechanics/molecular mechanics (QM/MM) computer simulations with transition path sampling and inertial likelihood maximization. We have identified a new solvent coordinate to complement the original atom-exchange coordinate used in the classic analysis by Chandrasekhar, Smith, and Jorgensen [J. Am. Chem. Soc. 107, 154 (1985)]. The new solvent coordinate quantifies instantaneous solvent-induced polarization relative to the equilibrium average charge density at each point along the reaction pathway. On the basis of likelihood scores and committor distributions, the new solvent coordinate improves upon the description of solvent dynamical effects relative to previously proposed solvent coordinates. However, it does not increase the transmission coefficient or the accuracy of a transition state theory rate calculation.

Diabat method for polymorph free energies: Extension to molecular crystals.

Lattice-switch Monte Carlo and the related diabat methods have emerged as efficient and accurate ways to compute free energy differences between polymorphs. In this work, we introduce a one-to-one mapping from the reference positions and displacements in one molecular crystal to the positions and displacements in another. Two features of the mapping facilitate lattice-switch Monte Carlo and related diabat methods for computing polymorph free energy differences. First, the mapping is unitary so that its Jacobian does not complicate the free energy calculations. Second, the mapping is easily implemented for molecular crystals of arbitrary complexity. We demonstrate the mapping by computing free energy differences between polymorphs of benzene and carbamazepine. Free energy calculations for thermodynamic cycles, each involving three independently computed polymorph free energy differences, all return to the starting free energy with a high degree of precision. The calculations thus provide a force field independent validation of the method and allow us to estimate the precision of the individual free energy differences.