Field-based Simulations Research Articles

Prediction of Thermochemical Properties of Long-Chain Alkanes Using Linear Regression: Application to Hydroisomerization.

Linear regression (LR) is used to predict thermochemical properties of alkanes at temperatures (0-1000) K to study chemical reaction equilibria inside zeolites. The thermochemical properties of C1 until C10 isomers reported by Scott are used as training data sets in the LR model which is used to predict these properties for alkanes longer than C10 isomers. Second-order groups are used as independent variables which account for the interactions between the neighboring groups of atoms. This model accurately predicts Gibbs free energies, enthalpies, Gibbs free energies of formation, and enthalpies of formation for alkanes which exceeds the chemical accuracy of 1 kcal/mol and outperforms the group contribution methods developed by Benson et al., Joback and Reid, and Constantinou and Gani. Predictions from our model are used to compute the reaction equilibrium distribution of hydroisomerization of C10 and C14 isomers in MTW-type zeolite. Calculation of reaction equilibrium distribution inside zeolites also requires Henry coefficients of the isomers which can be computed using classical force field-based molecular simulations using the RASPA2 software for which we created an automated workflow. The reaction equilibrium distribution for C10 isomers obtained using the LR model and the training data set for this model are in very good agreement. The tools developed in this study will enable the computational study of hydroisomerization of long-chain alkanes (>C10).

Photon Pathways and the Nonperturbative Scaling Law of High Harmonic Generation.

High harmonic generation (HHG) has become a core pillar of attosecond science. Traditionally described with field-based models, HHG can also be viewed in a parametric picture, which predicts all properties of the emitted photons, but not the nonperturbative efficiency of the process. Driving HHG with two noncollinear beams and deriving analytically the corresponding yield scaling laws for any intensity ratio, we herein reconcile the two interpretations, introducing a generalized photonic description of HHG. It is in full agreement with field-based simulations and experimental data, opening the route to smart engineering of HHG with multiple driving beams.

Combining particle and field-theoretic polymer models with multi-representation simulations.

Particle-based and field-theoretic simulations are both widely used methods to predict the properties of polymeric materials. In general, the advantages of each method are complementary. Field-theoretic simulations are preferred for polymers with high molecular weights and can provide direct access to chemical potentials and free energies, which makes them the method-of-choice for calculating phase diagrams. The trade-off is that field-theoretic simulations sacrifice the molecular details present in particle-based simulations, such as the configurations of individual molecules and their dynamics. In this work, we describe a new approach to conduct "multi-representation" simulations that efficiently map between particle-based and field-theoretic simulations. Our approach involves the construction of formally equivalent particle-based and field-based models, which are then simulated subject to the constraint that their spatial density profiles are equal. This constraint provides the ability to directly link particle-based and field-based simulations and enables calculations that can switch between one representation to the other. By switching between particle/field representations during a simulation, we demonstrate that our approach can leverage many of the advantages of each representation while avoiding their respective limitations. Although our method is illustrated in the context of complex sphere phases in linear diblock copolymers, we anticipate that it will be useful whenever free energies, rapid equilibration, molecular configurations, and dynamic information are all simultaneously desired.

Building polarization into protein-inhibitor binding dynamics in rational drug design for rheumatoid arthritis

Standard force field-based simulations to accomplish structure-based evaluations of lead molecules is a powerful tool. Combining protein fragmentation into tractable sub-systems with continuum solvation method is envisaged to enable quantum mechanics-based electronic structure calculations of macromolecules in their realistic environment. This along with incorporation of many-body polarization effect in molecular dynamics simulations may augment an accurate description of electrostatics of protein-inhibitor systems for effective drug design. Rheumatoid arthritis (RA) is a complex autoimmune disorder plagued by the ceiling effect of current targeted therapies, encouraging identification of new druggable targets and corresponding drug design to tackle the refractory form of disease. In this study, polarization-inclusive force field approach has been used to model protein solvation and ligand binding for ‘Mitogen-activated protein kinase’ (MAP3K8), a regulatory node of notable pharmacological relevance in RA synovial biology. For MAP3K8 inhibitors belonging to different scaffold series, the calculations illustrated differential electrostatic contribution to their relative binding affinities and successfully explained examples from available structure-activity relationship studies. Results from this study exemplified i) the advantage of this approach in reliably ranking inhibitors having close nanomolar range activities for the same target; and ii) its prospective application in lead molecule identification aiding drug discovery efforts in RA. Communicated by Ramaswamy H. Sarma

Accurate and Ultrafast Simulation of Molecular Recognition and Assembly on Metal Surfaces in Four Dimensions.

Understanding molecular interactions with metal surfaces in high reliability is critical for the development of catalysts, sensors, and therapeutics. Obtaining accurate experimental data for a wide range of surfaces remains a critical bottleneck and quantum-mechanical data remain speculative due to high uncertainties and limitations in scale. We report molecular dynamics simulations of adsorption energies and assembly of organic molecules on elemental metal surfaces using the INTERFACE force field (IFF). The force field-based simulations reach up to 8 times higher accuracy than density functional calculations at a million-fold faster speed, as well as more than 1 order of magnitude higher accuracy than other force fields relative to accurate measurements by single-crystal adsorption calorimetry. Uncertainties of prior computational methods are effectively reduced from on the order of 100% to less than 10% and validated by experimental data from multiple sources. Specifically, we describe the molecular interactions of benzene and naphthalene with even and defective platinum surfaces across a wide range of surface coverage in depth. We discuss molecular-scale influences on the heat of adsorption and clarify the definition of surface coverage. The methods can be applied to 18 metals to accurately predict binding and assembly of organic molecules, ligands, electrolytes, biological molecules, and gases without additional fit parameters.

Enhancement of formic acid production from carbon dioxide hydrogenation using metal-organic frameworks: Monte Carlo simulation study

Formic acid production from CO2 allows the reduction of carbon dioxide emissions while synthesizing a product with a wide range of applications. CO2 hydrogenation is challenging due to the cost of transition metal catalysts and the toxicity of the transition elements. In this work, the thermodynamic confinement effects of the metal–organic framework UiO-66 on the CO2 hydrogenation to formic acid were studied by force field-based molecular simulations. The confinement effects of UiO-66 and the metal–organic frameworks Cu-BTC, and IRMOF-1 were compared, to assess the impact of different pore size and metal centers on the production of HCOOH. Monte Carlo simulations in the grand-canonical ensemble were performed in the frameworks, using gas phase mole fractions of CO2, H2, and HCOOH at chemical equilibrium, obtained from Continuous Fractional Component Monte Carlo simulations in the Reaction Ensemble. The adsorption isobars of the components in metal–organic frameworks were computed at 298.15 – 800 K, 1 – 60 bar. The enhancement of HCOOH production due to preferential adsorption of HCOOH in metal–organic frameworks was calculated for all studied conditions. UiO-66, Cu-BTC, and IRMOF-1 affect CO2 hydrogenation reaction, shifting the thermodynamical equilibrium toward HCOOH formation. The prevailing factor is the type of metal center in the metal–organic framework. The confinement effect of Cu-BTC turns out to exceed the enhancement caused by UiO-66, and IRMOF-1. The resulting mole fraction of HCOOH increased by ca. 2000 times compared to the gas phase at 298.15 K, 60 bar. Cu-BTC can be considered as an alternative to improve the production of HCOOH due to elimination of the high-cost temperature elevation, cost reduction of downstream processing methods, and comparable final concentration of HCOOH to the reported concentrations of formate obtained using transition metal catalysts.

Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids.

Accurate empirical force fields of lipid molecules are a critical component of molecular dynamics simulation studies aimed at investigating properties of monolayers, bilayers, micelles, vesicles, and liposomes, as well as heterogeneous systems, such as protein-membrane complexes, bacterial cell walls, and more. While the majority of lipid force field-based simulations have been performed using pairwise-additive nonpolarizable models, advances have been made in the development of the polarizable force field based on the classical Drude oscillator model. In the present study, we undertake further optimization of the Drude lipid force field, termed Drude2023, including improved treatment of the phosphate and glycerol linker region of PC and PE headgroups, additional optimization of the alkene group in monounsaturated lipids, and inclusion of long-range Lennard-Jones interactions using the particle-mesh Ewald method. Initial optimization targeted quantum mechanical (QM) data on small model compounds representative of the linker region. Subsequent optimization targeted QM data on larger model compounds, experimental data, and dihedral potentials of mean force from the CHARMM36 additive lipid force field using a parameter reweighting protocol. The use of both experimental and QM target data during the reweighting protocol is shown to produce physically reasonable parameters that reproduce a collection of experimental observables. Target data for optimization included surface area/lipid for DPPC, DSPC, DMPC, and DLPC bilayers and nuclear magnetic resonance (NMR) order parameters for DPPC bilayers. Validation data include prediction of membrane thickness, scattering form factors, electrostatic potential profiles, compressibility moduli, surface area per lipid, water permeability, NMR T1 relaxation times, diffusion constants, and monolayer surface tensions for a variety of saturated and unsaturated lipid mono- and bilayers. Overall, the agreement with experimental data is quite good, though the results are less satisfactory for the NMR T1 relaxation times for carbons near the ester groups. Notable improvements compared to the additive C36 force field were obtained for membrane dipole potentials, lipid diffusion coefficients, and water permeability with the exception of monounsaturated lipid bilayers. It is anticipated that the optimized polarizable Drude2023 force field will help generate more accurate molecular simulations of pure bilayers and heterogeneous systems containing membranes, advancing our understanding of the role of electronic polarization in these systems.

Quantum Informed Machine-Learning Potentials for Molecular Dynamics Simulations of CO2's Chemisorption and Diffusion in Mg-MOF-74.

Among various porous solids for gas separation and purification, metal-organic frameworks (MOFs) are promising materials that potentially combine high CO2 uptake and CO2/N2 selectivity. So far, within the hundreds of thousands of MOF structures known today, it remains a challenge to computationally identify the best suited species. First principle-based simulations of CO2 adsorption in MOFs would provide the necessary accuracy; however, they are impractical due to the high computational cost. Classical force field-based simulations would be computationally feasible; however, they do not provide sufficient accuracy. Thus, the entropy contribution that requires both accurate force fields and sufficiently long computing time for sampling is difficult to obtain in simulations. Here, we report quantum-informed machine-learning force fields (QMLFFs) for atomistic simulations of CO2 in MOFs. We demonstrate that the method has a much higher computational efficiency (∼1000×) than the first-principle one while maintaining the quantum-level accuracy. As a proof of concept, we show that the QMLFF-based molecular dynamics simulations of CO2 in Mg-MOF-74 can predict the binding free energy landscape and the diffusion coefficient close to experimental values. The combination of machine learning and atomistic simulation helps achieve more accurate and efficient in silico evaluations of the chemisorption and diffusion of gas molecules in MOFs.

Measurements and simulations of light-activated matter.

Non-equilibrium active matter consumes energy to form self-organizing structures. For example, microtubules and motors assemble to form polar and nematic structures that facilitate cell division. Here, we investigate properties of varying aster geometries during and after contraction. Experiments, agent-based simulations, and field-based simulations provide insights into an in-vitro, light-activated system of microtubules (MTs) and kinesin motors. Using an agent-based simulation, we modeled the contraction process and computationally created radially symmetric aster structures; simulation results show asters contract faster at higher motor speeds and lower unbinding rates. Furthermore, the agent-based simulation captures 3D MT flux, potentially providing insight into field-based equations describing contraction. We aim to compare experimental contraction with a model from computational MT flux data, searching for mechanics that govern the formation and stabilization of aster structures.

Explainable Solvation Free Energy Prediction Combining Graph Neural Networks with Chemical Intuition.

The prediction of a molecule's solvation Gibbs free (ΔGsolv) energy in a given solvent is an important task which has traditionally been carried out via quantum chemical continuum methods or force field-based molecular simulations. Machine learning (ML) and graph neural networks in particular have emerged as powerful techniques for elucidating structure-property relationships. This work presents a graph neural network (GNN) for the prediction of ΔGsolv which, in addition to encoding typical atom and bond-level features, incorporates chemically intuitive, solvation-relevant parameters into the featurization process: semiempirical partial atomic charges and solvent dielectric constant. Solute-solvent interactions are included via an interaction map layer which can be visualized to examine solubility-enhancing or -decreasing interactions learnt by the model. On a test set of small organic molecules, our GNN predicts ΔGsolv in water and cyclohexane with an accuracy comparable to polarizable and ab initio generated force field methods [mean absolute error (MAE) = 0.4 and 0.2 kcal mol-1, respectively], without the need for any molecular simulation. For the FreeSolv data set of hydration free energies, the test MAE is 0.7 kcal mol-1. Interpretability and applicability of the model is highlighted through several examples including rationalizing the increased solubility of modified diaminoanthraquinones in organic solvents. The clear explanations afforded by our GNN allow for easy understanding of the model's predictions, giving the experimental chemist confidence in employing ML models toward more optimized synthetic routes.

CAT: A Compound Attachment Tool for the Construction of Composite Chemical Compounds.

The continuous improvement of computer architectures allows for the simulation of molecular systems of growing sizes. However, such calculations still require the input of initial structures, which are also becoming increasingly complex. In this work, we present CAT, a Compound Attachment Tool (source code available at https://github.com/nlesc-nano/CAT) and Python package for the automatic construction of composite chemical compounds, which supports the functionalization of organic, inorganic, and hybrid organic-inorganic materials. The CAT workflow consists in defining the anchoring sites on the reference material, usually a large molecular system denoted as a scaffold, and on the molecular species that are attached to it, i.e., the ligands. Usually, ligands are pre-optimized in a conformation biased toward more linear structures to minimize interligand(s) steric interactions, a bias that is important when multiple ligands are attached onto the scaffold. The resulting superstructure(s) are then stored in various formats that can be used afterward in quantum chemical calculations or classical force field-based simulations.

A simple molecular simulation strategy for rapid prediction of BTEX sorption on a surface-modified sorbent

Sorption of benzene, toluene, ethylbenzene, o-xylene, m-xylene, and p-xylene (BTEX) molecules on Fe3O4@SiO2@C18 sorbent allows the idea of entering the BTEX between the C18 (octadecyl) chains. This process resembles BTEX+C18 mixing. Given the dominant effect of non-bond interactions, it should be possible to predict the sorption behavior using force field-based molecular simulations. Experimental data, i.e. sorption efficiency of Fe3O4@SiO2@C18 towards BTEX in aqueous environment, were compared with results of two simulation strategies aimed at predicting sorption behavior using models significantly simplified compared to the real Fe3O4@SiO2@C18 structure. The first strategy involved molecular dynamics performed on models containing only the shell with C18 chains in water with BTEX. The second strategy involved miscibility calculations (based on modified Flory-Huggins theory) performed on models containing only pairs of molecules: C18+BTEX and H2O+BTEX. Results of both simulation strategies are in good agreement with experimental data, i.e. BTEX sorption on Fe3O4@SiO2@C18 can be studied using significantly simplified models. Given the speed of miscibility calculations and the simplicity of models used (pairs of molecules), the preparation of much larger models and time-consuming molecular dynamics simulations are not necessary. The sorption efficiency can be easily and quickly predicted by the miscibility calculations.

Condensation and growth of amorphous aluminosilicate nanoparticles via an aggregation process.

The precipitation of zeolite nanoparticles involves the initial formation of metastable precursors, such as amorphous entities, that crystallize through non-classical pathways. Here, using reactive force field-based simulations, we reveal how aluminosilicate oligomers grow concomitantly to the decondensation of silicate entities during the initial step of the reaction. Aluminate clusters first form in the solution, thus violating the Loewenstein rule in the first instant of the reaction, which is then followed by their connection with silicate oligomers at the terminal silanol groups before reorganization to finally diffuse within the silicate oligomers to form stable amorphous aluminosilicate nanoparticles that do obey the Loewenstein rule. Our results clearly indicate that aluminate does not serve as the nucleation center for the growth of aluminosilicates in a nucleation-like process but rather proceeds via an aggregation process. The coexistence of aluminosilicate oligomers and small silicate entities induces a phase separation that promotes the precipitation of zeolites with aging.

Efficient Parametrization of Force Field for the Quantitative Prediction of the Physical Properties of Ionic Liquid Electrolytes.

The prediction of transport properties of room-temperature ionic liquids from nonpolarizable force field-based simulations has long been a challenge. The uniform charge scaling method has been widely used to improve the agreement with the experiment by incorporating the polarizability and charge transfer effects in an effective manner. While this method improves the performance of the force fields, this prescription is ad hoc in character; further, a quantitative prediction is still not guaranteed. In such cases, the nonbonded interaction parameters too need to be refined, which requires significant effort. In this work, we propose a three-step semiautomated refinement procedure based on (1) atomic site charges obtained from quantum calculations of the bulk condensed phase; (2) quenched Monte Carlo optimizer to shortlist suitable force field candidates, which are then tested using pilot simulations; and (3) manual refinement to further improve the accuracy of the force field. The strategy is designed in a sequential manner with each step improving the accuracy over the previous step, allowing the users to invest the effort commensurate with the desired accuracy of the refined force field. The refinement procedure is applied on N,N-diethyl-N-methyl-N-(2-methoxyethyl)ammonium bis(trifluoromethanesulfonyl)imide (DEME-TFSI), a front-runner as an electrolyte for electric double-layer capacitors and single-molecule-based devices. The transferability of the refined force field is tested on N,N-dimethyl-N-ethyl-N-methoxyethoxyethylammonium bis(trifluoromethanesulfonyl)imide (N112,2O2O1-TFSI). The refined force field is found to be better at predicting both structural and transport properties compared to the uniform charge scaling procedure, which showed a discrepancy in the X-ray structure factor. The refined force field showed quantitative agreement with structural (density and X-ray structure factor) and transport properties-diffusion coefficients, ionic conductivity, and shear viscosity over a wide temperature range, building a case for the wide adoption of the procedure.

Molecularly Informed Field Theories from Bottom-up Coarse-Graining.

Polymer formulations possessing mesostructures or phase coexistence are challenging to simulate using atomistic particle-explicit approaches due to the disparate time and length scales, while the predictive capability of field-based simulations is hampered by the need to specify interactions at a coarser scale (e.g., χ-parameters). To overcome the weaknesses of both, we introduce a bottom-up coarse-graining methodology that leverages all-atom molecular dynamics to molecularly inform coarser field-theoretic models. Specifically, we use relative-entropy coarse-graining to parametrize particle models that are directly and analytically transformable into statistical field theories. We demonstrate the predictive capability of this approach by reproducing experimental aqueous poly(ethylene oxide) (PEO) cloud-point curves with no parameters fit to experimental data. This synergistic approach to multiscale polymer simulations opens the door to de novo exploration of phase behavior across a wide variety of polymer solutions and melt formulations.

Effect of Water Content on Thermodynamic Properties of Compressed Hydrogen.

Force field-based molecular simulations were used to calculate thermal expansivities, heat capacities, and Joule–Thomson coefficients of binary (standard) hydrogen–water mixtures for temperatures between 366.15 and 423.15 K and pressures between 50 and 1000 bar. The mole fraction of water in saturated hydrogen–water mixtures in the gas phase ranges from 0.004 to 0.138. The same properties were calculated for pure hydrogen at 323.15 K and pressures between 100 and 1000 bar. Simulations were performed using the TIP3P and a modified TIP4P force field for water and the Marx, Vrabec, Cracknell, Buch, and Hirschfelder force fields for hydrogen. The vapor–liquid equilibria of hydrogen–water mixtures were calculated along the melting line of ice Ih, corresponding to temperatures between 264.21 and 272.4 K, using the TIP3P force field for water and the Marx force field for hydrogen. In this temperature range, the solubilities and the chemical potentials of hydrogen and water were obtained. Based on the computed solubility data of hydrogen in water, the freezing-point depression of water was computed ranging from 264.21 to 272.4 K. The modified TIP4P and Marx force fields were used to improve the solubility calculations of hydrogen–water mixtures reported in our previous study [RahbariA.;J. Chem. Eng. Data2019, 64, 4103−4115] for temperatures between 323 and 423 K and pressures ranging from 100 to 1000 bar. The chemical potentials of ice Ih were calculated as a function of pressure between 100 and 1000 bar, along the melting line for temperatures between 264.21 and 272.4 K, using the IAPWS equation of state for ice Ih. We show that at low pressures, the presence of water has a large effect on the thermodynamic properties of compressed hydrogen. Our conclusions may have consequences for the energetics of a hydrogen refueling station using electrochemical hydrogen compressors.

PySpawn: Software for Nonadiabatic Quantum Molecular Dynamics.

The ab initio multiple spawning (AIMS) method enables nonadiabatic quantum molecular dynamics simulations in an arbitrary number of dimensions, with potential energy surfaces provided by electronic structure calculations performed on-the-fly. However, the intricacy of the AIMS algorithm complicates software development, deployment on modern shared computer resources, and postsimulation data analysis. PySpawn is a nonadiabatic molecular dynamics software package that addresses these issues. The program is designed to be easily interfaced with electronic structure software, and an interface to the TeraChem software package is described here. PySpawn introduces a task-based reorganization of the AIMS algorithm, allowing fine-grained restart capability and setting the stage for efficient parallelization in a future release. PySpawn includes a user-friendly and interactive Python analysis module that will enable novice users to painlessly adopt AIMS. As a demonstration of PySpawn's simulation capability and analysis module, we report complete active space self-consistent field-based AIMS simulations of the 1,2-dithienyl-1,2-dicyanoethene molecule, a promising molecular photoswitch.

Two-dimensional diamine-linked covalent organic frameworks for CO2/N2 capture and separation: theoretical modeling and simulations.

Two-dimensional covalent organic frameworks (2D-COFs) with diamine-based linkers have been designed and investigated for CO2/N2 gaseous mixture adsorption and separation via a systematic theoretical study by combining density functional theory (DFT) calculations and force field-based molecular dynamics (MD) simulations. We explored the adsorption sites and adsorption energies of CO2/N2 on 2D-COFs. The gas uptake capacity, adsorption isotherms, permeability, and selectivity were simulated based on an improved formulation of force fields for mixture separation in post-combustion conditions. This theoretical approach provided atomistic understanding and quantitative description of intermolecular interactions governing the physisorption dynamics of the considered systems. The results suggest that 2D-COFs investigated in this study are competitive with other 2D materials for carbon capture and separation and can be considered as alternative molecular sieving materials offering efficient and rapid separation and adsorption of different molecules.

Defects, Diffusion, and Dopants in Li2Ti6O13: Atomistic Simulation Study.

In this study, force field-based simulations are employed to examine the defects in Li-ion diffusion pathways together with activation energies and a solution of dopants in Li2Ti6O13. The lowest defect energy process is found to be the Li Frenkel (0.66 eV/defect), inferring that this defect process is most likely to occur. This study further identifies that cation exchange (Li–Ti) disorder is the second lowest defect energy process. Long-range diffusion of Li-ion is observed in the bc-plane with activation energy of 0.25 eV, inferring that Li ions move fast in this material. The most promising trivalent dopant at the Ti site is Co3+, which would create more Li interstitials in the lattice required for high capacity. The favorable isovalent dopant is the Ge4+ at the Ti site, which may alter the mechanical property of this material. The electronic structures of the favorable dopants are analyzed using density functional theory (DFT) calculations.

Phase Behavior of Grafted Polymer Nanocomposites from Field-Based Simulations

There are limited theoretically predicted phase diagrams for polymer nanocomposites (PNCs) because conventional modeling techniques are largely unable to predict the macroscale phase behavior of PN...