Building Interpretable Machine Learning Models to Identify Chemometric Trends in Seabirds of the North Pacific Ocean.

Marine environmental monitoring efforts often rely on the bioaccumulation of persistent anthropogenic contaminants in organisms to create a spatiotemporal record of the ecosystem. Intercorrelation results from the origin, uptake, and transport of these contaminants throughout the ecosystem and may be affected by organism-specific processes such as biotransformation. Here, we explore trends that machine learning tools reveal about a large, recently released environmental chemistry data set of common anthropogenic pollutants measured in the eggs of five seabird species from the North Pacific Ocean. We modeled these data with a variety of machine learning approaches and found models that could accurately determine a range of taxonomic and spatiotemporal trends. We illustrate a general workflow and set of analysis tools that can be used to identify interpretable models which perform nearly as well as state-of-the-art "black boxes." For example, we found shallow decision trees that could resolve genus with greater than 96% accuracy using as few as two analytes and a k-nearest neighbor classifier that could resolve species differences with more than 94% accuracy using only five analytes. The benefits of interpretability outweighed the marginally improved accuracy of more complex models. This demonstrates how machine learning may be used to discover rational, quantitative trends in these systems.

Derivable genetic programming for two-dimensional colloidal materials.

We describe a method for deriving surface functionalization patterns for colloidal systems that can induce self-assembly into any chosen periodic symmetry at a planar interface. The result is a sequence of letters, s ∈ {A,T,C,G}, or a gene, that describes the perimeter of the colloidal object and programs its self-assembly. This represents a genome that is finite and can be exhaustively enumerated. These genes derive from symmetry, which may be topologically represented by two-dimensional parabolic orbifolds; since these orbifolds are surfaces that may be derived from first principles, this represents an ab initio route to colloid functionality. The genes are human readable and can be employed to easily design colloidal units. We employ a biological (genetic) analogy to demonstrate this and illustrate their connection to the designs of Maurits Cornelis (M. C.) Escher.

Symmetry-derived structure directing agents for two-dimensional crystals of arbitrary colloids.

We derive properties of self-assembling rings which can template the organization of an arbitrary colloid into any periodic symmetry in two Euclidean dimensions. By viewing this as a tiling problem, we illustrate how the shape and chemical patterning of these rings are derivable, and are explicitly reflected by the symmetry group's orbifold symbol. We performed molecular dynamics simulations to observe their self-assembly and found 5 different characteristics which could be easily rationalized on the basis of this symbol. These include systems which undergo chiral phase separation, are addressably complex, exhibit self-limiting growth into clusters, form ordered "rods" in only one-dimension akin to a smectic phase, and those from symmetry groups which are pluripotent and allow one to select rings which exhibit different behaviors. We discuss how the curvature of the ring's edges plays an integral role in achieving correct self-assembly, and illustrate how to obtain these shapes. This provides a method for patterning colloidal systems at interfaces without explicitly programming this information onto the colloid itself.

Grand canonical inverse design of multicomponent colloidal crystals.

Inverse design methods are powerful computational approaches for creating colloidal systems which self-assemble into a target morphology by reverse engineering the Hamiltonian of the system. Despite this, these optimization procedures tend to yield Hamiltonians which are too complex to be experimentally realized. An alternative route to complex structures involves the use of several different components, however, conventional inverse design methods do not explicitly account for the possibility of phase separation into compositionally distinct structures. Here, we present an inverse design scheme for multicomponent colloidal systems by combining active learning with a method to directly compute their ground state phase diagrams. This explicitly accounts for phase separation and can locate stable regions of Hamiltonian parameter space which grid-based surveys are prone to miss. Using this we design low-density, binary structures with Lennard-Jones-like pairwise interactions that are simpler than in the single component case and potentially realizable in an experimental setting. This reinforces the concept that ground states of simple, multicomponent systems might be rich with previously unappreciated diversity, enabling the assembly of non-trivial structures with only few simple components instead of a single complex one.

Dynamic arrest of adhesive hard rod dispersions.

The phenomenon of dynamic arrest, more commonly referred to as gel and glass formation, originates as particle motion slows significantly. Current understanding of gels and glasses stems primarily from dispersions of spherical particles, but much less is known about how particle shape affects dynamic arrest transitions. To better understand the effects of particle shape anisotropy on gel and glass formation, we systematically measure the rheology, particle dynamics, and static microstructure of thermoreversible colloidal dispersions of adhesive hard rods (AHR). First, the dynamic arrest transitions are mapped as a function of temperature T, aspect ratio L/D≈ 3 to 7, and volume fraction φ≈ 0.1 to 0.5. The critical gel temperature Tgel and glass volume fraction φg are determined from the particle dynamics and rheology. Second, an effective orientation-averaged, short-range attraction between rods is quantified from small-angle scattering measurements and characterized by a reduced temperature τ. Similar τ is found at low rod concentrations, indicating that rod gelation occurs at similar effective attraction strength independent of L/D. Monte Carlo simulations reveal a similar convergence in τ when rods cluster and percolate with an average bond coordination number 〈nc〉≈ 2.4, supporting the link between physical gelation and rigidity percolation. Lastly, AHR results are mapped onto a dimensionless state diagram to compare with previous predictions of attraction-driven gels, repulsion-driven glasses, and liquid crystal phases.

Symmetry-Based Crystal Structure Enumeration in Two Dimensions.

The accurate prediction of stable crystalline phases is a long-standing problem encountered in the study of conventional atomic and molecular solids as well as soft materials. One possible solution involves enumerating a reasonable set of candidate structures and then screening them to identify the one(s) with the lowest (free) energy. Candidate structures in this set can also serve as starting points for other routines, such as genetic algorithms, which search via optimization. Here, we present a framework for crystal structure enumeration of two-dimensional systems that utilizes a combination of symmetry- and stoichiometry-imposed constraints to compute valid configurations of particles that tile Euclidean space. With mild assumptions, this produces a computationally tractable total number of proposed candidates, enabling multicomponent systems to be screened by direct enumeration of possible crystalline ground states. The python code that enables these calculations is available at https://github.com/usnistgov/PACCS.

Extrapolation and interpolation strategies for efficiently estimating structural observables as a function of temperature and density.

Thermodynamic extrapolation has previously been used to predict arbitrary structural observables in molecular simulations at temperatures (or relative chemical potentials in open-system mixtures) different from those at which the simulation was performed. This greatly reduces the computational cost in mapping out phase and structural transitions. In this work, we explore the limitations and accuracy of thermodynamic extrapolation applied to water, where qualitative shifts from anomalous to simple-fluid-like behavior are manifested through shifts in the liquid structure that occur as a function of both temperature and density. We present formulas for extrapolating in volume for canonical ensembles and demonstrate that linear extrapolations of water's structural properties are only accurate over a limited density range. On the other hand, linear extrapolation in temperature can be accurate across the entire liquid state. We contrast these extrapolations with classical perturbation theory techniques, which are more conservative and slowly converging. Indeed, we show that such behavior is expected by demonstrating exact relationships between extrapolation of free energies and well-known techniques to predict free energy differences. An ideal gas in an external field is also studied to more clearly explain these results for a toy system with fully analytical solutions. We also present a recursive interpolation strategy for predicting arbitrary structural properties of molecular fluids over a predefined range of state conditions, demonstrating its success in mapping qualitative shifts in water structure with density.

On the virial coefficients of confined fluids: Analytic expressions for the thermodynamic properties of hard particles with attractions in slit and cylindrical pores to second order.

We present analytic expressions for the second virial coefficient of hard particles in slit and cylindrical pores with arbitrary square-well fluid-fluid and fluid-solid interactions, number of fluid components, and pore sizes. We derive multiple methods to predict the isotherms in confined materials based on the virial expansion and find that the thermodynamic and structural properties calculated from the virial expansion to second order in the dilute limit display excellent qualitative agreement with previous simulation results.

Assembly of three-dimensional binary superlattices from multi-flavored particles.

Binary superlattices constructed from nano- or micron-sized colloidal particles have a wide variety of applications, including the design of advanced materials. Self-assembly of such crystals from their constituent colloids can be achieved in practice by, among other means, the functionalization of colloid surfaces with single-stranded DNA sequences. However, when driven by DNA, this assembly is traditionally premised on the pairwise interaction between a single DNA sequence and its complement, and often relies on particle size asymmetry to entropically control the crystalline arrangement of its constituents. The recently proposed "multi-flavoring" motif for DNA functionalization, wherein multiple distinct strands of DNA are grafted in different ratios to different colloids, can be used to experimentally realize a binary mixture in which all pairwise interactions are independently controllable. In this work, we use various computational methods, including molecular dynamics and Wang-Landau Monte Carlo simulations, to study a multi-flavored binary system of micron-sized DNA-functionalized particles modeled implicitly by Fermi-Jagla pairwise interactions. We show how self-assembly of such systems can be controlled in a purely enthalpic manner, and by tuning only the interactions between like particles, demonstrate assembly into various morphologies. Although polymorphism is present over a wide range of pairwise interaction strengths, we show that careful selection of interactions can lead to the generation of pure compositionally ordered crystals. Additionally, we show how the crystal composition changes with the like-pair interaction strengths, and how the solution stoichiometry affects the assembled structures.

Molecular simulation of capillary phase transitions in flexible porous materials.

We used flat-histogram sampling Monte Carlo to study capillary phase transitions in deformable adsorbent materials. Specifically, we considered a pure adsorbate fluid below its bulk critical temperature within a slit pore of variable pore width. The instantaneous pore width is dictated by a number of factors, such as adsorbate loading, reservoir pressure, fluid-wall interaction, and bare adsorbent properties. In the slit pores studied here, the bare adsorbent free energy was assumed to be biparabolic, consisting of two preferential pore configurations, namely, the narrow pore and the large pore configurations. Four distinct phases could be found in the adsorption isotherms. We found a low-pressure phase transition, driven primarily by capillary condensation/evaporation and accompanied by adsorbent deformation in response. The deformation can be a relatively small contraction/expansion as seen in elastic materials, or a large-scale structural transformation of the adsorbent. We also found a high-pressure transition driven by excluded volume effects, which tends to expand the material and thus results in a large-scale structural transformation of the adsorbent. The adsorption isotherms and osmotic free energies can be rationalized by considering the relative free energy differences between the basins of the bare adsorbent free energy.

Predicting structural properties of fluids by thermodynamic extrapolation.

We describe a methodology for extrapolating the structural properties of multicomponent fluids from one thermodynamic state to another. These properties generally include features of a system that may be computed from an individual configuration such as radial distribution functions, cluster size distributions, or a polymer's radius of gyration. This approach is based on the principle of using fluctuations in a system's extensive thermodynamic variables, such as energy, to construct an appropriate Taylor series expansion for these structural properties in terms of intensive conjugate variables, such as temperature. Thus, one may extrapolate these properties from one state to another when the series is truncated to some finite order. We demonstrate this extrapolation for simple and coarse-grained fluids in both the canonical and grand canonical ensembles, in terms of both temperatures and the chemical potentials of different components. The results show that this method is able to reasonably approximate structural properties of such fluids over a broad range of conditions. Consequently, this methodology may be employed to increase the computational efficiency of molecular simulations used to measure the structural properties of certain fluid systems, especially those used in high-throughput or data-driven investigations.

Assembly of multi-flavored two-dimensional colloidal crystals.

We systematically investigate the assembly of binary multi-flavored colloidal mixtures in two dimensions. In these mixtures all pairwise interactions between species may be tuned independently. This introduces an additional degree of freedom over more traditional binary mixtures with fixed mixing rules, which is anticipated to open new avenues for directed self-assembly. At present, colloidal self-assembly into non-trivial lattices tends to require either high pressures for isotropically interacting particles, or the introduction of directionally anisotropic interactions. Here we demonstrate tunable assembly into a plethora of structures which requires neither of these conditions. We develop a minimal model that defines a three-dimensional phase space containing one dimension for each pairwise interaction, then employ various computational techniques to map out regions of this phase space in which the system self-assembles into these different morphologies. We then present a mean-field model that is capable of reproducing these results for size-symmetric mixtures, which reveals how to target different structures by tuning pairwise interactions, solution stoichiometry, or both. Concerning particle size asymmetry, we find that domains in this model's phase space, corresponding to different morphologies, tend to undergo a continuous "rotation" whose magnitude is proportional to the size asymmetry. Such continuity enables one to estimate the relative stability of different lattices for arbitrary size asymmetries. Owing to its simplicity and accuracy, we expect this model to serve as a valuable design tool for engineering binary colloidal crystals from multi-flavored components.

Tuning flexibility to control selectivity in soft porous crystals.

We use flat-histogram Monte Carlo simulations to study how changing the flexibility of soft porous crystals (SPCs) affects their selective adsorption of a binary, size-asymmetric supercritical fluid. Specifically, we consider mesoporous SPCs which have multiple minima in their free energy profiles as a function of pore size such that they are capable of exhibiting polymorphism between a narrow and large pore phase. While specific fluid-pore interactions determine the shape of both pores' selectivity curve as a function of adsorbate pressure, an individual pore tends to selectively adsorb a species based on the size of the adsorbate molecule relative to itself, thereby shifting the pore's selectivity curve relative to its polymorph. By controlling the flexibility of a SPC, the relative thermodynamic stability of the two pore phases may be varied, thereby changing the overall selectivity of the SPC during adsorbate loading. We investigate this for two classes of SPCs: one representative of "gate-opening" materials and another of "breathing" materials. For gate-opening materials, this control is much more salient than in breathing ones. However, for the latter, we illustrate how to tune the free energy profile to create materials which breathe multiple times during adsorption/desorption.

Controlling relative polymorph stability in soft porous crystals with a barostat.

We use Monte Carlo simulations to investigate the thermodynamic behavior of soft porous crystal (SPC) adsorbents under the influence of an external barostat. We consider SPCs that naturally exhibit polymorphism between crystal forms of two distinct pore sizes. In the absence of barostatting, these crystals may be naturally divided into two categories depending on their response to stress applied by the adsorbate fluid: those which macroscopically deform and change the volume of their unit cell ("breathing") and those which instead undergo internal rearrangements that change the adsorbate-accessible volume without modifying the unit cell volume ("gate-opening"). When breathing SPCs have a constant external pressure applied, in addition to the thermodynamic pressure of the adsorbate fluid, we find that the free energy difference between the crystal polymorphs is shifted by a constant amount over the entire course of adsorption. Thus, their relative stability may be easily controlled by the barostat. However, when the crystal is held at a fixed overall pressure, changes to the relative stability of the polymorphs tend to be more complex. We demonstrate a thermodynamic analogy between breathing SPCs held at a fixed pressure and macroscopically rigid gate-opening ones which explains this behavior. Furthermore, we illustrate how this implies that external mechanical forces may be employed to tune the effective free energy profile of an empty SPC, which may open new avenues to engineer the thermodynamic properties of these polymorphic adsorbents, such as selectivity.

Temperature extrapolation of multicomponent grand canonical free energy landscapes.

We derive a method for extrapolating the grand canonical free energy landscape of a multicomponent fluid system from one temperature to another. Previously, we introduced this statistical mechanical framework for the case where kinetic energy contributions to the classical partition function were neglected for simplicity [N. A. Mahynski et al., J. Chem. Phys. 146, 074101 (2017)]. Here, we generalize the derivation to admit these contributions in order to explicitly illustrate the differences that result. Specifically, we show how factoring out kinetic energy effects a priori, in order to consider only the configurational partition function, leads to simpler mathematical expressions that tend to produce more accurate extrapolations than when these effects are included. We demonstrate this by comparing and contrasting these two approaches for the simple cases of an ideal gas and a non-ideal, square-well fluid.

Multivariable extrapolation of grand canonical free energy landscapes.

We derive an approach for extrapolating the free energy landscape of multicomponent systems in the grand canonical ensemble, obtained from flat-histogram Monte Carlo simulations, from one set of temperature and chemical potentials to another. This is accomplished by expanding the landscape in a Taylor series at each value of the order parameter which defines its macrostate phase space. The coefficients in each Taylor polynomial are known exactly from fluctuation formulas, which may be computed by measuring the appropriate moments of extensive variables that fluctuate in this ensemble. Here we derive the expressions necessary to define these coefficients up to arbitrary order. In principle, this enables a single flat-histogram simulation to provide complete thermodynamic information over a broad range of temperatures and chemical potentials. Using this, we also show how to combine a small number of simulations, each performed at different conditions, in a thermodynamically consistent fashion to accurately compute properties at arbitrary temperatures and chemical potentials. This method may significantly increase the computational efficiency of biased grand canonical Monte Carlo simulations, especially for multicomponent mixtures. Although approximate, this approach is amenable to high-throughput and data-intensive investigations where it is preferable to have a large quantity of reasonably accurate simulation data, rather than a smaller amount with a higher accuracy.

Predicting low-temperature free energy landscapes with flat-histogram Monte Carlo methods.

We present a method for predicting the free energy landscape of fluids at low temperatures from flat-histogram grand canonical Monte Carlo simulations performed at higher ones. We illustrate our approach for both pure and multicomponent systems using two different sampling methods as a demonstration. This allows us to predict the thermodynamic behavior of systems which undergo both first order and continuous phase transitions upon cooling using simulations performed only at higher temperatures. After surveying a variety of different systems, we identify a range of temperature differences over which the extrapolation of high temperature simulations tends to quantitatively predict the thermodynamic properties of fluids at lower ones. Beyond this range, extrapolation still provides a reasonably well-informed estimate of the free energy landscape; this prediction then requires less computational effort to refine with an additional simulation at the desired temperature than reconstruction of the surface without any initial estimate. In either case, this method significantly increases the computational efficiency of these flat-histogram methods when investigating thermodynamic properties of fluids over a wide range of temperatures. For example, we demonstrate how a binary fluid phase diagram may be quantitatively predicted for many temperatures using only information obtained from a single supercritical state.

Communication: Predicting virial coefficients and alchemical transformations by extrapolating Mayer-sampling Monte Carlo simulations.

Virial coefficients are predicted over a large range of both temperatures and model parameter values (i.e., alchemical transformation) from an individual Mayer-sampling Monte Carlo simulation by statistical mechanical extrapolation with minimal increase in computational cost. With this extrapolation method, a Mayer-sampling Monte Carlo simulation of the SPC/E (extended simple point charge) water model quantitatively predicted the second virial coefficient as a continuous function spanning over four orders of magnitude in value and over three orders of magnitude in temperature with less than a 2% deviation. In addition, the same simulation predicted the second virial coefficient if the site charges were scaled by a constant factor, from an increase of 40% down to zero charge. This method is also shown to perform well for the third virial coefficient and the exponential parameter for a Lennard-Jones fluid.

Multicomponent adsorption in mesoporous flexible materials with flat-histogram Monte Carlo methods.

We demonstrate an extensible flat-histogram Monte Carlo simulation methodology for studying the adsorption of multicomponent fluids in flexible porous solids. This methodology allows us to easily obtain the complete free energy landscape for the confined fluid-solid system in equilibrium with a bulk fluid of any arbitrary composition. We use this approach to study the adsorption of a prototypical coarse-grained binary fluid in "Hookean" solids, where the free energy of the solid may be described as a simple spring. However, our approach is fully extensible to solids with arbitrarily complex free energy profiles. We demonstrate that by tuning the fluid-solid interaction ranges, the inhomogeneous fluid structure inside the pore can give rise to enhanced selective capture of a larger species through cooperative adsorption with a smaller one. The maximum enhancement in selectivity is observed at low to intermediate pressures and is especially pronounced when the larger species is very dilute in the bulk. This suggest a mechanism by which the selective capture of a minor component from a bulk fluid may be enhanced.

Relative stability of the FCC and HCP polymorphs with interacting polymers.

Recent work [Mahynski et al., Nat. Commun., 2014, 5, 4472] has demonstrated that the addition of long linear homopolymers thermodynamically biases crystallizing hard-sphere colloids to produce the hexagonal close-packed (HCP) polymorph over the closely related face-centered cubic (FCC) structure when the polymers and colloids are purely repulsive. In this report, we investigate the effects of thermal interactions on each crystal polymorph to explore the possibility of stabilizing the FCC crystal structure over the HCP. We find that the HCP polymorph remains at least as stable as its FCC counterpart across the entire range of interactions we explored, where interactions were quantified by the reduced second virial coefficient, -1.50 < B < 1.01. This metric conveniently characterizes the crossover from entropically to energetically dominated systems at B ≈ 0. While the HCP relies on its octahedral void arrangement for enhanced stability when B > 0, its tetrahedral voids produce a similar effect when B < 0 (i.e. when energetics dominate). Starting from this, we derive a mean-field expression for the free energy of an infinitely-dilute polymer adsorbed in the crystal phase for nonzero B. Our results reveal that co-solute biasing of a single polymorph can still be observed in experimentally realizable scenarios when the colloids and polymers have attractive interactions, and provide a possible explanation for the experimental finding that pure FCC crystals are elusive in these binary mixtures.