MoSDeF-GOMC: Python Software for the Creation of Scientific Workflows for the Monte Carlo Simulation Engine GOMC.

MoSDeF-GOMC is a python interface for the Monte Carlo software GOMC to the Molecular Simulation Design Framework (MoSDeF) ecosystem. MoSDeF-GOMC automates the process of generating initial coordinates, assigning force field parameters, and writing coordinate (PDB), connectivity (PSF), force field parameter, and simulation control files. The software lowers entry barriers for novice users while allowing advanced users to create complex workflows that encapsulate simulation setup, execution, and data analysis in a single script. All relevant simulation parameters are encoded within the workflow, ensuring reproducible simulations. MoSDeF-GOMC's capabilities are illustrated through a number of examples, including prediction of the adsorption isotherm for CO2 in IRMOF-1, free energies of hydration for neon and radon over a broad temperature range, and the vapor-liquid coexistence curve of a four-component surrogate for the jet fuel S-8. The MoSDeF-GOMC software is available on GitHub at https://github.com/GOMC-WSU/MoSDeF-GOMC.

Transferable Force Fields from Experimental Scattering Data with Machine Learning Assisted Structure Refinement.

Deriving transferable pair potentials from experimental neutron and X-ray scattering measurements has been a longstanding challenge in condensed matter physics. State-of-the-art scattering analysis techniques estimate real-space microstructure from reciprocal-space total scattering data by refining pair potentials to obtain agreement between simulated and experimental results. Prior attempts to apply these potentials with molecular simulations have revealed inaccurate predictions of thermodynamic fluid properties. In this Letter, a machine learning assisted structure-inversion method applied to neutron scattering patterns of the noble gases (Ne, Ar, Kr, and Xe) is shown to recover transferable pair potentials that accurately reproduce both microstructure and vapor-liquid equilibria from the triple to critical point. Therefore, it is concluded that a single neutron scattering measurement is sufficient to predict macroscopic thermodynamic properties over a wide range of states and provide novel insight into local atomic forces in dense monatomic systems.

Self-Assembly and Biogenesis of the Cellular Membrane are Dictated by Membrane Stretch and Composition.

The cell plasma membrane is a highly dynamic organelle governing a wide range of cellular activities including ion transport, secretion, cell division, growth, and development. The fundamental process involved in the addition of new membranes to pre-existing plasma membranes, however, is unclear. Here, we report, using biophysical, morphological, biochemical, and molecular dynamic simulations, the selective incorporation of proteins and lipids from the cytosol into the cell plasma membrane dictated by membrane stretch and composition. Stretching of the cell membrane as a consequence of volume increase following incubation in a hypotonic solution and results in the incorporation of cytosolic proteins and lipids into the existing plasma membrane. Molecular dynamic simulations further confirm that increased membrane stretch results in the rapid insertion of lipids into the existing plasma membrane. Similarly, depletion of cholesterol from the cell plasma membrane selectively alters the incorporation of lipids into the membrane.

Molecular exchange Monte Carlo: A generalized method for identity exchanges in grand canonical Monte Carlo simulations.

A generalized identity exchange algorithm is presented for Monte Carlo simulations in the grand canonical ensemble. The algorithm, referred to as molecular exchange Monte Carlo, may be applied to multicomponent systems of arbitrary molecular topology and provides significant enhancements in the sampling of phase space over a wide range of compositions and temperatures. Three different approaches are presented for the insertion of large molecules, and the pros and cons of each method are discussed. The performance of the algorithms is highlighted through grand canonical Monte Carlo histogram-reweighting simulations performed on a number of systems, which include methane+n-alkanes, butane+perfluorobutane, water+impurity, and 2,2,4-trimethylpentane+neopentane. Relative acceptance efficiencies for molecule transfers of up to 400 times that of standard configurational-bias Monte Carlo are obtained.

Optimized Mie potentials for phase equilibria: Application to noble gases and their mixtures with n-alkanes.

Transferrable force fields, based on n-6 Mie potentials, are presented for noble gases. By tuning the repulsive exponent, ni, it is possible to simultaneously reproduce experimental saturated liquid densities and vapor pressures with high accuracy, from the normal boiling point to the critical point. Vapor-liquid coexistence curves for pure fluids are calculated using histogram reweighting Monte Carlo simulations in the grand canonical ensemble. For all noble gases, saturated liquid densities and vapor pressures are reproduced to within 1% and 4% of experiment, respectively. Radial distribution functions, extracted from NVT and NPT Monte Carlo simulations, are in similarly excellent agreement with experimental data. The transferability of the optimized force fields is assessed through calculations of binary mixture vapor-liquid equilibria. These mixtures include argon + krypton, krypton + xenon, methane + krypton, methane + xenon, krypton + ethane, and xenon + ethane. For all mixtures, excellent agreement with experiment is achieved without the introduction of any binary interaction parameters or multi-body interactions.

Biomolecular simulations with the transferable potentials for phase equilibria: extension to phospholipids.

The Transferable Potentials for Phase Equilibria (TraPPE) is extended to zwitterionic and charged lipids including phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and phosphatidylglycerol (PG). The performance of the force field is validated through isothermal-isobaric ensemble (NPT) molecular dynamics simulations of hydrated lipid bilayers performed with the aforementioned head groups combined with saturated and unsaturated alkyl tails containing 12-18 carbon atoms. The effects of water model and sodium ion parameters on the performance of the lipid force field are determined. The predictions of the TraPPE force field for the area per lipid, bilayer thickness, and volume per lipid are within 1-5% of experimental values. Key structural properties of the bilayer, such as order parameter splitting in the sn-2 chain and X-ray form factors, are found to be in close agreement with experimental data.

Prediction of 1-octanol-water and air-water partition coefficients for nitro-aromatic compounds from molecular dynamics simulations.

United-atom force fields, based on the Transferable Potentials for Phase Equilibria (TraPPE), are developed for twelve nitro-aromatic compounds, which include 2,4-dinitrotoluene (2,4-DNT), 2,6-dinitrotoluene (2,6-DNT), 3-nitrotoluene (3-NT), 4-nitrotoluene (4-NT), 1,3-dinitrobenzene (1,3-DNB), 1,4-dinitrobenzene (1,4-DNB), 2,4-dinitroanisole (DNAN), 1,3,5-trinitrobenzene (TNB), 2,4,6-trinitrotoluene (TNT), 2-nitroanisole (2-NAN), 4-nitroanisole (4-NAN) and n-methyl-p-nitroaniline (MNA). 1-Octanol-water and air-water partition coefficients are predicted for the optimized TraPPE-UA force field with adaptive biasing force molecular dynamics simulations, and compared to available experimental data. Log Kow values are predicted with an average absolute deviation of 0.2 log units, while Henry's law constants are predicted to with an average absolute deviation of 0.5 log units. Two additional models are presented for energetic materials with five membered rings for which no experimental data are available in the open literature: 3,5-dinitropyrazole (DNP) and 3-nitro-1,2,4-triazole-5-one (NTO). Investigation of the local microstructure around each solute reveals that 1-octanol is able to form hydrogen bonded chains around the solute, while little organized microstructure was observed around the solutes in water.

Direct calculation of 1-octanol-water partition coefficients from adaptive biasing force molecular dynamics simulations.

The 1-octanol-water partition coefficient log K(ow) of a solute is a key parameter used in the prediction of a wide variety of complex phenomena such as drug availability and bioaccumulation potential of trace contaminants. In this work, adaptive biasing force molecular dynamics simulations are used to determine absolute free energies of hydration, solvation, and 1-octanol-water partition coefficients for n-alkanes from methane to octane. Two approaches are evaluated; the direct transfer of the solute from 1-octanol to water phase, and separate transfers of the solute from the water or 1-octanol phase to vacuum, with both methods yielding statistically indistinguishable results. Calculations performed with the TIP4P and SPC∕E water models and the TraPPE united-atom force field for n-alkanes show that the choice of water model has a negligible effect on predicted free energies of transfer and partition coefficients for n-alkanes. A comparison of calculations using wet and dry octanol phases shows that the predictions for log K(ow) using wet octanol are 0.2-0.4 log units lower than for dry octanol, although this is within the statistical uncertainty of the calculation.

Computational prediction of ionic liquid 1-octanol/water partition coefficients.

Wet 1-octanol/water partition coefficients (log K(ow)) predicted for imidazolium-based ionic liquids using adaptive bias force-molecular dynamics (ABF-MD) simulations lie in excellent agreement with experimental values. These encouraging results suggest prospects for this computational tool in the a priori prediction of log K(ow) values of ionic liquids broadly with possible screening implications as well (e.g., prediction of CO(2)-philic ionic liquids).

Monte Carlo predictions of phase equilibria and structure for dimethyl ether + sulfur dioxide and dimethyl ether + carbon dioxide.

A new force field for dimethyl ether (DME) based on the Lennard-Jones (LJ) 12-6 plus point charge functional form is presented in this work. This force field reproduces experimental saturated liquid and vapor densities, vapor pressures, heats of vaporization, and critical properties to within the statistical uncertainty of the combined experimental and simulation measurements for temperatures between the normal boiling and critical point. Critical parameters and normal boiling point are predicted to within 0.1% of experiment. This force field is used in grand canonical histogram reweighting Monte Carlo simulations to predict the pressure composition diagrams for the binary mixtures DME + SO(2) at 363.15 K and DME + CO(2) at 335.15 and 308.15 K. For the DME + SO(2) mixture, simulation is able to qualitatively reproduce the minimum pressure azeotropy observed experimentally for this mixture, but quantitative errors exist, suggesting that multibody effects may be important in this system. For the DME + CO(2) mixture, simulation is able to predict the pressure-composition behavior within 1% of experimental data. Simulations in the isobaric-isothermal ensemble are used to determine the microstructure of DME + SO(2) and DME + CO(2) mixtures. The DME + SO(2) shows weak pairing between DME and SO(2) molecules, while no specific pairing or aggregation is observed for mixtures of DME + CO(2).

Development of an optimized intermolecular potential for sulfur dioxide.

A new force field for sulfur dioxide, capable of predicting accurately the vapor-liquid equilibria, critical properties, vapor pressure, and heats of vaporization is presented. The new force field reproduces the saturated liquid densities, vapor pressures and heats of vaporization to within 0.5, 2, and 2% of experiment, respectively. The predicted critical properties and the normal boiling point are in excellent agreement with experimental results. Pair distribution functions are calculated for the S-S, S-O, and O-O interactions are in close agreement with neutron and X-ray scattering experiments. In addition to the new force field, similar calculations are performed for four SO(2) intermolecular potentials proposed by Sokolic et al. (Sokolic, F.; Guissani, Y. and Guillot, B. J. Phys. Chem. 1985, 89, 3023], which show that these models work reasonably well near the state point where they were originally parametrized, but large errors in the predicted coexistence properties are displayed at higher and lower temperatures. Comparison of the radial distribution functions show the local structure is only weakly affected by the different force field parameters.

Ca(2+) bridging of apposed phospholipid bilayers.

In an effort to provide insight into the mechanism of Ca(2+)-induced fusion of lipid vesicles, molecular dynamics simulations in the isobaric-isothermal ensemble are used to investigate interactions of Ca(2+) with apposed lipid bilayers in close proximity. Simulations reveal the formation of a Ca(2+)-phospholipid "anhydrous complex" between apposed bilayers, whereas similar calculations performed with Na(+) display only complexation between neighboring lipids within the same bilayer. The binding of Ca(2+) to apposed phospholipids brings large regions of the bilayers into close contact (<4 Å), displacing water from phospholipid head groups in the process and creating regions of local dehydration. Dehydration of the apposed bilayers leads to ordering of the phospholipid tails, which is partially disrupted by the presence of Ca(2+)-phospholipid bridges.

Mie potentials for phase equilibria calculations: application to alkanes and perfluoroalkanes.

Transferable united-atom force fields, based on n - 6 Lennard-Jones potentials, are presented for normal alkanes and perfluorocarbons. It is shown that by varying the repulsive exponent the range of the potential can be altered, leading to improved predictions of vapor pressures while also reproducing saturated liquid densities to high accuracy. Histogram-reweighting Monte Carlo simulations in the grand canonical ensemble are used to determine the vapor liquid coexistence curves, vapor pressures, heats of vaporization, and critical points for normal alkanes methane through tetradecane, and perfluorocarbons perfluoromethane through perfluorooctane. For all molecules studied, saturated liquid densities are reproduced to within 1% of experiment. Vapor pressures for normal alkanes and perfluorocarbons were predicted to within 3% and 6% of experiment, respectively. Calculations performed for binary mixture vapor-liquid equilibria for propane + pentane show excellent agreement with experiment, while slight deviations are observed for the ethane + perfluoroethane mixture.

Extension of the transferable potentials for phase equilibria force field to dimethylmethyl phosphonate, sarin, and soman.

The transferable potentials for phase equilibria force field is extended to dimethylmethylphosphonate (DMMP), sarin, and soman by introducing a new interaction site representing the phosphorus atom. Parameters for the phosphorus atom are optimized to reproduce the liquid densities at 303 and 373 K and the normal boiling point of DMMP. Calculations for sarin and soman are performed in predictive mode, without further parameter optimization. Vapor-liquid coexistence curves, critical properties, vapor pressures and heats of vaporization are predicted over a wide range of temperatures with histogram reweighting Monte Carlo simulations in the grand canonical ensemble. Excellent agreement with experiment is achieved for all compounds, with unsigned errors of less than 1% for vapor pressures and normal boiling points and under 5% for heats of vaporization and liquid densities at ambient conditions.

All-atom force field for the prediction of vapor-liquid equilibria and interfacial properties of HFA134a.

A new all-atom force field capable of accurately predicting the bulk and interfacial properties of 1,1,1,2-tetrafluoroethane (HFA134a) is reported. Parameterization of several force fields with different initial charge configurations from ab initio calculations was performed using the histogram reweighting method and Monte Carlo simulations in the grand canonical ensemble. The 12-6 Lennard-Jones well depth and diameter for the different HFA134a models were determined by fitting the simulation results to pure-component vapor-equilibrium data. Initial screening of the force fields was achieved by comparing the calculated and experimental bulk properties. The surface tension of pure HFA134a served as an additional screening property to help discriminate an optimum model. The proposed model reproduces the experimental saturated liquid and vapor densities, and the vapor pressure for HFA134a within average errors of 0.7%, 4.4%, and 3.1%, respectively. Critical density, temperature, vapor pressure, normal boiling point, and heat of vaporization at 298 K are also in good agreement with experimental data with errors of 0.2%, 0.1%, 6.2%, 0%, 2.2%, respectively. The calculated surface tension is found to be within the experimental range of 7.7-8.1 mN.m(-1). The dipole moment of the different models was found to significantly affect the prediction of the vapor pressure and surface tension. The ability of the HFA134a models in predicting the interfacial tension against water is also discussed. The results presented here are relevant in the development of technologies where the more environmentally friendly HFA134a is utilized as a substitute to the ozone depleting chlorofluorocarbon propellants.

Ca2+-dimethylphosphate complex formation: providing insight into Ca2+-mediated local dehydration and membrane fusion in cells.

Earlier studies using X-ray diffraction, light scattering, photon correlation spectroscopy, and atomic force microscopy, strongly suggest that SNARE-induced membrane fusion in cells proceeds as a result of calcium bridging opposing bilayers. The bridging of phospholipid heads groups in the opposing bilayers by calcium leads to the release of water from hydrated Ca(2+) ions as well as the loosely coordinated water at PO-lipid head groups. Local dehydration of phospholipid head groups and the calcium, bridging opposing bilayers, then leads to destabilization of the lipid bilayers and membrane fusion. This hypothesis was tested in the current study by atomistic molecular dynamic simulations in the isobaric-isothermal ensemble using hydrated dimethylphosphate anions (DMP(-)) and calcium cations. Results from the study demonstrate, formation of DMP-Ca(2+) complexes and the consequent removal of water, supporting the hypothesis. Our study further demonstrates that as a result of Ca(2+)-DMP self-assembly, the distance between anionic oxygens between the two DMP molecules is reduced to 2.92A, which is in close agreement with the 2.8A SNARE-induced apposition established between opposing bilayers, reported earlier from X-ray diffraction measurements.

Transferable potentials for phase equilibria. 8. United-atom description for thiols, sulfides, disulfides, and thiophene.

An extension of the transferable potentials for phase equilibria united-atom (TraPPE-UA) force field to thiol, sulfide, and disulfide functionalities and thiophene is presented. In the TraPPE-UA force field, nonbonded interactions are governed by a Lennard-Jones plus fixed point charge functional form. Partial charges are determined through a CHELPG analysis of electrostatic potential energy surfaces derived from ab initio calculations at the HF/6-31g+(d,p) level. The Lennard-Jones well depth and size parameters for four new interaction sites, S (thiols), S(sulfides), S(disulfides), and S(thiophene), were determined by fitting simulation data to pure-component vapor-equilibrium data for methanethiol, dimethyl sulfide, dimethyl disulfide, and thiophene, respectively. Configurational-bias Monte Carlo simulations in the grand canonical ensemble combined with histogram-reweighting methods were used to calculate the vapor-liquid coexistence curves for methanethiol, ethanethiol, 2-methyl-1-propanethiol, 2-methyl-2-propanethiol, 2-butanethiol, pentanethiol, octanethiol, dimethyl sulfide, diethyl sulfide, ethylmethyl sulfide, dimethyl disulfide, diethyl disulfide, and thiophene. Excellent agreement with experiment is achieved, with unsigned errors of less than 1% for saturated liquid densities and less than 3% for critical temperatures. The normal boiling points were predicted to within 1% of experiment in most cases, although for certain molecules (pentanethiol) deviations as large as 5% were found. Additional calculations were performed to determine the pressure-composition behavior of ethanethiol+n-butane at 373.15 K and the temperature-composition behavior of 1-propanethiol+n-hexane at 1.01 bar. In each case, a good reproduction of experimental vapor-liquid equilibrium separation factors is achieved; both of the coexistence curves are somewhat shifted because of overprediction of the pure-component vapor pressures.

Molecular modeling of phase behavior and microstructure of acetone-chloroform-methanol binary mixtures.

Force fields based on a Lennard-Jones (LJ) 12-6 plus point charge functional form are developed for acetone and chloroform specifically to reproduce the minimum pressure azeotropy found experimentally in this system. Point charges are determined from a CHELPG population analysis performed on an acetone-chloroform dimer. The required electrostatic surface for this dimer is determined from ab initio calculations performed with MP2 theory and the 6-31g++(3df,3pd) basis set. LJ parameters are then optimized such that the liquid-vapor coexistence curve, critical parameters, and vapor pressures are well reproduced by simulation. Histogram-reweighting Monte Carlo simulations in the grand canonical ensemble are used to determine the phase diagrams for the binary mixtures acetone-chloroform, acetone-methanol, and chloroform-methanol. The force fields developed in this work reproduce the minimum pressure azeotrope in the acetone-chloroform mixture found in experiment. The predicted azeotropic composition of x(CHCl3) = 0.77 is in fair agreement with the experimental value of x(CHCl3)expt = 0.64. The new force fields were also found to provide improved predictions of the pressure-composition behavior of acetone-methanol and chloroform-methanol when compared to other force fields commonly used for vapor-liquid equilibria calculations. NPT simulations were conducted at 300 K and 1 bar for equimolar mixtures of acetone-chloroform, acetone-methanol, and methanol-chloroform. Analysis of the microstructure reveals significant hydrogen bonding occurring between acetone and chloroform. Limited interspecies hydrogen bonding was found in the acetone-methanol or chloroform-methanol mixtures.

Effect of partial charge parametrization on the fluid phase behavior of hydrogen sulfide.

The effect of partial charge parametrization on the fluid phase behavior of hydrogen sulfide is investigated with grand canonical histogram reweighting Monte Carlo simulations. Four potential models, based on a Lennard-Jones + point charge functional form, are developed. It is shown that Lennard-Jones parameters can be tuned such that partial charges for the sulfur atom in the range -0.40 < q(s) < -0.252 will lead to an accurate reproduction of experimental vapor-liquid equilibria. Each of the parameter sets developed in this work are used to predict the pressure composition behavior H2S-n-pentane at 377.6 K. While the mixture calculation provides a means of reducing the number of candidate parameter sets, multiple parameter sets were found to yield an excellent reproduction of both the pure component and mixture phase behavior.