Monte Carlo molecular simulations with FEASST version 0.25.1.

FEASST is an open-source Monte Carlo software package for particle-based simulations. This software, which was released in 2017, has been used to study phase equilibrium, self-assembly, aggregation or gelation in biological materials, colloids, polymers, ionic liquids, and adsorption in porous networks. We highlight some of the unique features available in FEASST, such as flat-histogram grand canonical ensemble, Gibbs ensemble, and Mayer-sampling simulations with support for anisotropic models and parallelization with flat-histogram and prefetching. We also discuss how the challenges of supporting a variety of Monte Carlo algorithms were overcome by an object-oriented design. This also allows others to extend classes, which improves software interoperability, as inspired by LAMMPS classes and user packages. This article describes version 0.25.1 with benchmarks, compilation instructions, and introductory tutorials for running, restarting, and testing simulations, user guidelines, software design strategies, alternative interfaces, and the test-driven development strategy.

Anisotropic coarse-grain Monte Carlo simulations of lysozyme, lactoferrin, and NISTmAb by precomputing atomistic models.

We develop a multiscale coarse-grain model of the NIST Monoclonal Antibody Reference Material 8671 (NISTmAb) to enable systematic computational investigations of high-concentration physical instabilities such as phase separation, clustering, and aggregation. Our multiscale coarse-graining strategy captures atomic-resolution interactions with a computational approach that is orders of magnitude more efficient than atomistic models, assuming the biomolecule can be decomposed into one or more rigid bodies with known, fixed structures. This method reduces interactions between tens of thousands of atoms to a single anisotropic interaction site. The anisotropic interaction between unique pairs of rigid bodies is precomputed over a discrete set of relative orientations and stored, allowing interactions between arbitrarily oriented rigid bodies to be interpolated from the precomputed table during coarse-grained Monte Carlo simulations. We present this approach for lysozyme and lactoferrin as a single rigid body and for the NISTmAb as three rigid bodies bound by a flexible hinge with an implicit solvent model. This coarse-graining strategy predicts experimentally measured radius of gyration and second osmotic virial coefficient data, enabling routine Monte Carlo simulation of medically relevant concentrations of interacting proteins while retaining atomistic detail. All methodologies used in this work are available in the open-source software Free Energy and Advanced Sampling Simulation Toolkit.

Theory and Monte Carlo simulation of the ideal gas with shell particles in the canonical, isothermal-isobaric, grand canonical, and Gibbs ensembles.

Theories of small systems play an important role in the fundamental understanding of finite size effects in statistical mechanics, as well as the validation of molecular simulation results as no computer can simulate fluids in the thermodynamic limit. Previously, a shell particle was included in the isothermal-isobaric ensemble in order to resolve an ambiguity in the resulting partition function. The shell particle removed either redundant volume states or redundant translational degrees of freedom of the system and yielded quantitative differences from traditional simulations in this ensemble. In this work, we investigate the effect of including a shell particle in the canonical, grand canonical, and Gibbs ensembles. For systems comprised of a pure component ideal gas, analytical expressions for various thermodynamic properties are obtained. We also derive the Metropolis Monte Carlo simulation acceptance criteria for these ensembles with shell particles, and the results of the simulations of an ideal gas are in excellent agreement with the theoretical predictions. The system size dependence of various important ensemble averages is also analyzed.

Flow Activation Energy of High-Concentration Monoclonal Antibody Solutions and Protein-Protein Interactions Influenced by NaCl and Sucrose.

The solution viscosity and protein-protein interactions (PPIs) as a function of temperature (4-40 °C) were measured at a series of protein concentrations for a monoclonal antibody (mAb) with different formulation conditions, which include NaCl and sucrose. The flow activation energy (Eη) was extracted from the temperature dependence of solution viscosity using the Arrhenius equation. PPIs were quantified via the protein diffusion interaction parameter (kD) measured by dynamic light scattering, together with the osmotic second virial coefficient and the structure factor obtained through small-angle X-ray scattering. Both viscosity and PPIs were found to vary with the formulation conditions. Adding NaCl introduces an attractive interaction but leads to a significant reduction in the viscosity. However, adding sucrose enhances an overall repulsive effect and leads to a slight decrease in viscosity. Thus, the averaged (attractive or repulsive) PPI information is not a good indicator of viscosity at high protein concentrations for the mAb studied here. Instead, a correlation based on the temperature dependence of viscosity (i.e., Eη) and the temperature sensitivity in PPIs was observed for this specific mAb. When kD is more sensitive to the temperature variation, it corresponds to a larger value of Eη and thus a higher viscosity in concentrated protein solutions. When kD is less sensitive to temperature change, it corresponds to a smaller value of Eη and thus a lower viscosity at high protein concentrations. Rather than the absolute value of PPIs at a given temperature, our results show that the temperature sensitivity of PPIs may be a more useful metric for predicting issues with high viscosity of concentrated solutions. In addition, we also demonstrate that caution is required in choosing a proper protein concentration range to extract kD. In some excipient conditions studied here, the appropriate protein concentration range needs to be less than 4 mg/mL, remarkably lower than the typical concentration range used in the literature.

Flat-Histogram Monte Carlo Simulation of Water Adsorption in Metal-Organic Frameworks.

Molecular simulations of water adsorption in porous materials often converge slowly due to sampling bottlenecks that follow from hydrogen bonding and, in many cases, the formation of water clusters. These effects may be exacerbated in metal-organic framework (MOF) adsorbents, due to the presence of pore spaces (cages) that promote the formation of discrete-size clusters and hydrophobic effects (if present), among other reasons. In Grand Canonical Monte Carlo (MC) simulations, these sampling challenges are typically manifested by low MC acceptance ratios, a tendency for the simulation to become stuck in a particular loading state (i.e., macrostates), and the persistence of specific clusters for long periods of the simulation. We present simulation strategies to address these sampling challenges, by applying flat-histogram MC (FHMC) methods and specialized MC move types to simulations of water adsorption. FHMC, in both Transition-matrix and Wang-Landau forms, drives the simulation to sample relevant macrostates by incorporating weights that are self-consistently adjusted throughout the simulation and generate the macrostate probability distribution (MPD). Specialized MC moves, based on aggregation-volume bias and configurational bias methods, separately address low acceptance ratios for basic MC trial moves and specifically target water molecules in clusters; in turn, the specialized MC moves improve the efficiency of generating new configurations which is ultimately reflected in improved statistics collected by FHMC. The combined strategies are applied to study the adsorption of water in CuBTC and ZIF-8 at 300 K, through examination of the MPD and the adsorption isotherm generated by histogram reweighting. A key result is the appearance of nontrivial oscillations in the MPD, which we show to be associated with water clusters in the adsorption system. Additionally, we show that the probabilities of certain clusters become similar in value near the boundaries of the isotherm hysteresis loop, indicating a strong connection between cluster formation/destruction and the thermodynamic limits of stability. For a hydrophobic MOF, the FHMC results show that the phase transition from low density to high density is suppressed to water pressure far above the bulk-fluid saturation pressure; this is consistent with results presented elsewhere. We also compare our FHMC simulation isotherm to one measured by a different technique but with ostensibly the same molecular interactions and comment on observed differences and the need for follow-up work. The simulation strategies presented here can be applied to the simulation of water in other MOFs using heuristic guidelines laid out in our text, which should facilitate the more consistent and efficient simulation of water adsorption in porous materials in future applications.

Programming Interfacial Porosity and Symmetry with Escherized Colloids.

We simultaneously designed the porosity and plane symmetry of self-assembling colloidal films by using isohedral tiles to determine the location and shape of enthalpically interacting surface patches on motifs being functionalized. The symmetries of both the tile and motif determine the plane symmetry group of the final assembly. Previous work has either ignored symmetry considerations altogether or accounted for only the tile's properties, applicable only when the motif is asymmetric; this approach provides a complete account and enables the design of symmetric colloids using this tile-based approach, which are often more practical to manufacture. We present the methodology, based on the type of the tile, and provide computational tools that enable the automatic classification of all tiles for a given motif and the optimization of the tile to fit the motif, sometimes referred to as "Escherization".

Role of Domain-Domain Interactions on the Self-Association and Physical Stability of Monoclonal Antibodies: Effect of pH and Salt.

Monoclonal antibodies (mAbs) make up a major class of biotherapeutics with a wide range of clinical applications. Their physical stability can be affected by various environmental factors. For instance, an acidic pH can be encountered during different stages of the mAb manufacturing process, including purification and storage. Therefore, understanding the behavior of flexible mAb molecules in acidic solution environments will benefit the development of stable mAb products. This study used small-angle X-ray scattering (SAXS) and complementary biophysical characterization techniques to investigate the conformational flexibility and protein-protein interactions (PPI) of a model mAb molecule under near-neutral and acidic conditions. The study also characterized the interactions between Fab and Fc fragments under the same buffer conditions to identify domain-domain interactions. The results suggest that solution pH significantly influences mAb flexibility and thus could help mAbs remain physically stable by maximizing local electrostatic repulsions when mAbs become crowded in solution. Under acidic buffer conditions, both Fab and Fc contribute to the repulsive PPI observed among the full mAb at a low ionic strength. However, as ionic strength increases, hydrophobic interactions lead to the self-association of Fc fragments and, subsequently, could affect the aggregation state of the mAb.

Leveraging uncertainty estimates and derivative information in Gaussian process regression for efficient collection and use of molecular simulation data.

We introduce Gaussian Process Regression (GPR) as an enhanced method of thermodynamic extrapolation and interpolation. The heteroscedastic GPR models that we introduce automatically weight provided information by its estimated uncertainty, allowing for the incorporation of highly uncertain, high-order derivative information. By the linearity of the derivative operator, GPR models naturally handle derivative information and, with appropriate likelihood models that incorporate heterogeneous uncertainties, are able to identify estimates of functions for which the provided observations and derivatives are inconsistent due to the sampling bias that is common in molecular simulations. Since we utilize kernels that form complete bases on the function space to be learned, the estimated uncertainty in the model takes into account that of the functional form itself, in contrast to polynomial interpolation, which explicitly assumes the functional form to be fixed. We apply GPR models to a variety of data sources and assess various active learning strategies, identifying when specific options will be most useful. Our active-learning data collection based on GPR models incorporating derivative information is finally applied to tracing vapor-liquid equilibrium for a single-component Lennard-Jones fluid, which we show represents a powerful generalization to previous extrapolation strategies and Gibbs-Duhem integration. A suite of tools implementing these methods is provided at https://github.com/usnistgov/thermo-extrap.

Efficiency Comparison of Single- and Multiple-Macrostate Grand Canonical Ensemble Transition-Matrix Monte Carlo Simulations.

Recent interest in parallelizing flat-histogram transition-matrix Monte Carlo simulations in the grand canonical ensemble, due to its demonstrated effectiveness in studying phase behavior, self-assembly and adsorption, has led to the most extreme case of single-macrostate simulations, where each macrostate is simulated independently with ghost particle insertions and deletions. Despite their use in several studies, no efficiency comparisons of these single-macrostate simulations have been made with multiple-macrostate simulations. We show that multiple-macrostate simulations are up to 3 orders of magnitude more efficient than single-macrostate simulations, which demonstrates the remarkable efficiency of flat-histogram biased insertions and deletions, even with low acceptance probabilities. Efficiency comparisons were made for supercritical fluids and vapor-liquid equilibrium of bulk Lennard-Jones and a three-site water model, self-assembling patchy trimer particles and adsorption of a Lennard-Jones fluid confined in a purely repulsive porous network, using the open source simulation toolkit FEASST. By directly comparing with a variety of Monte Carlo trial move sets, this efficiency loss in single-macrostate simulations is attributed to three related reasons. First, ghost particle insertions and deletions in single-macrostate simulations incur the same computational expense as grand canonical ensemble trials in multiple-macrostate simulations, yet ghost trials do not reap the sampling benefit from propagating the Markov chain to a new microstate. Second, single-macrostate simulations lack macrostate change trials that are biased by the self-consistently converging relative macrostate probability, which is a major component of flat histogram simulations. Third, limiting a Markov chain to a single macrostate reduces sampling possibilities. Existing parallelization methods for multiple-macrostate flat-histogram simulations are shown to be more efficient than parallel single-macrostate simulations by approximately an order of magnitude or more in all systems investigated.

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.

Systematic control of collective variables learned from variational autoencoders.

Variational autoencoders (VAEs) are rapidly gaining popularity within molecular simulation for discovering low-dimensional, or latent, representations, which are critical for both analyzing and accelerating simulations. However, it remains unclear how the information a VAE learns is connected to its probabilistic structure and, in turn, its loss function. Previous studies have focused on feature engineering, ad hoc modifications to loss functions, or adjustment of the prior to enforce desirable latent space properties. By applying effectively arbitrarily flexible priors via normalizing flows, we focus instead on how adjusting the structure of the decoding model impacts the learned latent coordinate. We systematically adjust the power and flexibility of the decoding distribution, observing that this has a significant impact on the structure of the latent space as measured by a suite of metrics developed in this work. By also varying weights on separate terms within each VAE loss function, we show that the level of detail encoded can be further tuned. This provides practical guidance for utilizing VAEs to extract varying resolutions of low-dimensional information from molecular dynamics and Monte Carlo simulations.

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.

Temperature Extrapolation of Henry's Law Constants and the Isosteric Heat of Adsorption.

Computational screening of adsorbent materials often uses the Henry's law constant (KH) (at a particular temperature) as a first discriminator metric due to its relative ease of calculation. The isosteric heat of adsorption in the limit of zero pressure (qst∞) is often calculated along with the Henry's law constant, and both properties are informative metrics of adsorbent material performance at low-pressure conditions. In this article, we introduce a method for extrapolating KH as a function of temperature, using series-expansion coefficients that are easily computed at the same time as KH itself; the extrapolation function also yields qst∞. The extrapolation is highly accurate over a wide range of temperatures when the basis temperature is sufficiently high, for a wide range of adsorbent materials and adsorbate gases. Various results suggest that the extrapolation is accurate when the extrapolation range in inverse-temperature space is limited to |ß - ß0 | < 0.5 mol/kJ. Application of the extrapolation to a large set of materials is shown to be successful provided that KH is not extremely large and/or the extrapolation coefficients converge satisfactorily. The extrapolation is also able to predict qst∞ for a system that shows an unusually large temperature dependence. The work provides a robust method for predicting KH and qst∞ over a wide range of industrially relevant temperatures with minimal effort beyond that necessary to compute those properties at a single temperature, which facilitates the addition of practical operating (or processing) conditions to computational screening exercises.

Learning Efficient, Collective Monte Carlo Moves with Variational Autoencoders.

Discovering meaningful collective variables for enhancing sampling, via applied biasing potentials or tailored MC move sets, remains a major challenge within molecular simulation. While recent studies identifying collective variables with variational autoencoders (VAEs) have focused on the encoding and latent space discovered by a VAE, the impact of the decoding and its ability to act as a generative model remains unexplored. We demonstrate how VAEs may be used to learn (on-the-fly and with minimal human intervention) highly efficient, collective Monte Carlo moves that accelerate sampling along the learned collective variable. In contrast to many machine learning-based efforts to bias sampling and generate novel configurations, our methods result in exact sampling in the ensemble of interest and do not require reweighting. In fact, we show that the acceptance rates of our moves approach unity for a perfect VAE model. While this is never observed in practice, VAE-based Monte Carlo moves still enhance sampling of new configurations. We demonstrate, however, that the form of the encoding and decoding distributions, in particular the extent to which the decoder reflects the underlying physics, greatly impacts the performance of the trained VAE.

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.

Predicting the Mixing Behavior of Aqueous Solutions Using a Machine Learning Framework.

The most direct approach to determining if two aqueous solutions will phase-separate upon mixing is to exhaustively screen them in a pair-wise fashion. This is a time-consuming process that involves preparation of numerous stock solutions, precise transfer of highly concentrated and often viscous solutions, exhaustive agitation to ensure thorough mixing, and time-sensitive monitoring to observe the presence of emulsion characteristics indicative of phase separation. Here, we examined the pair-wise mixing behavior of 68 water-soluble compounds by observing the formation of microscopic phase boundaries and droplets of 2278 unique 2-component solutions. A series of machine learning classifiers (artificial neural network, random forest, k-nearest neighbors, and support vector classifier) were then trained on physicochemical property data associated with the 68 compounds and used to predict their miscibility upon mixing. Miscibility predictions were then compared to the experimental observations. The random forest classifier was the most successful classifier of those tested, displaying an average receiver operator characteristic area under the curve of 0.74. The random forest classifier was validated by removing either one or two compounds from the input data, training the classifier on the remaining data and then predicting the miscibility of solutions involving the removed compound(s) using the classifier. The accuracy, specificity, and sensitivity of the random forest classifier were 0.74, 0.80, and 0.51, respectively, when one of the two compounds to be examined was not represented in the training data. When asked to predict the miscibility of two compounds, neither of which were represented in the training data, the accuracy, specificity, and sensitivity values for the random forest classifier were 0.70, 0.82 and 0.29, respectively. Thus, there is potential for this machine learning approach to improve the design of screening experiments to accelerate the discovery of aqueous two-phase systems for numerous scientific and industrial applications.

On the virial expansion of model adsorptive systems.

We investigate the thermodynamic properties of various super-critical model adsorptive systems with different fluid-solid attractive strengths using the confined-density virial expansion, with coefficients calculated using the Mayer-sampling Monte Carlo method up to fifth order. We find that the virial expansion converges for adsorptive systems over a density range corresponding approximately to the film-formation regime. Beyond this regime, higher order effects become increasingly important. The virial expansion of the density profile is also investigated. It is determined that this expansion gives insight into the structure associated with adsorption. We also find that weakly attractive systems have a more negative second virial coefficient than strongly attractive systems. This runs counter to the usual interpretation of bulk fluid virial coefficients. This is due to the infinite-dilution limit being very different for adsorbed fluids compared to bulk fluids.

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.

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.