Simulation Of Large Systems Research Articles

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials.

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Role of Distributed Inverter-Based-Resources in Bulk Grid Primary Frequency Response Through HELICS Based SMTD Co-Simulation

In the recent years, there has been an increase in distributed inverter-based resources (D-IBRs) in the system that is impacting the grid performance. Dynamic transmission and distribution (T&D) co-simulation is an effective method to evaluate the impact of D-IBRs on the system response, however, implementation is a challenge. This article innovatively achieves accurate T&D model coordination with inclusion of dynamic models of power electronic controls. It also demonstrates that D-IBRs operating with grid-forming (GFM) controls with frequency support capability can aid in primary frequency regulation in a low-inertia power system. The article discusses the fundamentals of dynamic T&D co-simulation to analyze the primary frequency response on the IEEE 39-Bus transmission system coupled to IEEE 123-Bus distribution system with D-IBRs with GFM controls. To demonstrate scalability, a large system (286383 T&D buses) is co-simulated: IEEE 39-bus transmission system with all load buses coupled to a total of 2328 parallel IEEE 123-Bus distribution feeders containing a mix of D-IBRs is considered. The article also discusses the use of parallel computation used for distribution system solutions enabled in the scalable multi-timescale T&D co-simulation framework along with computation time comparisons between the large and small system simulations with GFM inverter dynamics.

Species-selective nanoreactor molecular dynamics simulations based on linear-scaling tight-binding quantum chemical calculations.

Here, extensions to quantum chemical nanoreactor molecular dynamics simulations for discovering complex reactive events are presented. The species-selective algorithm, where the nanoreactor effectively works for the selected desired reactants, was introduced to the original scheme. Moreover, for efficient simulations of large model systems with the modified approach, the divide-and-conquer linear-scaling density functional tight-binding method was exploited. Two illustrative applications of the polymerization of propylene and cyclopropane mixtures and the aggregation of sodium chloride from aqueous solutions indicate that species-selective quantum chemical nanoreactor molecular dynamics is a promising method to accelerate the sampling of multicomponent chemical processes proceeding under relatively mild conditions.

Exact non-adiabatic coupling vectors for the time-dependent density functional based tight-binding method.

We report on non-adiabatic coupling vectors between electronic excited states for the time-dependent-density functional theory based tight-binding (TD-DFTB) method. The implementation includes orbital relaxation effects that have been previously neglected and covers also the case of range-separated exchange-correlation functionals. Benchmark calculations with respect to first principles TD-DFT highlight the large dependence of non-adiabatic couplings on the functional. Closer investigations of the topology around a conical intersection between excited states show that TD-DFTB delivers near-exact values of the Berry phase, which paves the way for consistent non-adiabatic molecular dynamics simulations for large systems.

TBMaLT, a flexible toolkit for combining tight-binding and machine learning.

Tight-binding approaches, especially the Density Functional Tight-Binding (DFTB) and the extended tight-binding schemes, allow for efficient quantum mechanical simulations of large systems and long-time scales. They are derived from ab initio density functional theory using pragmatic approximations and some empirical terms, ensuring a fine balance between speed and accuracy. Their accuracy can be improved by tuning the empirical parameters using machine learning techniques, especially when information about the local environment of the atoms is incorporated. As the significant quantum mechanical contributions are still provided by the tight-binding models, and only short-ranged corrections are fitted, the learning procedure is typically shorter and more transferable as it were with predicting the quantum mechanical properties directly with machine learning without an underlying physically motivated model. As a further advantage, derived quantum mechanical quantities can be calculated based on the tight-binding model without the need for additional learning. We have developed the open-source framework-Tight-Binding Machine Learning Toolkit-which allows the easy implementation of such combined approaches. The toolkit currently contains layers for the DFTB method and an interface to the GFN1-xTB Hamiltonian, but due to its modular structure and its well-defined interfaces, additional atom-based schemes can be implemented easily. We are discussing the general structure of the framework, some essential implementation details, and several proof-of-concept applications demonstrating the perspectives of the combined methods and the functionality of the toolkit.

Regularized by Physics: Graph Neural Network Parametrized Potentials for the Description of Intermolecular Interactions.

Simulations of molecular systems using electronic structure methods are still not feasible for many systems of biological importance. As a result, empirical methods such as force fields (FF) have become an established tool for the simulation of large and complex molecular systems. The parametrization of FF is, however, time-consuming and has traditionally been based on experimental data. Recent years have therefore seen increasing efforts to automatize FF parametrization or to replace FF with machine-learning (ML) based potentials. Here, we propose an alternative strategy to parametrize FF, which makes use of ML and gradient-descent based optimization while retaining a functional form founded in physics. Using a predefined functional form is shown to enable interpretability, robustness, and efficient simulations of large systems over long time scales. To demonstrate the strength of the proposed method, a fixed-charge and a polarizable model are trained on ab initio potential-energy surfaces. Given only information about the constituting elements, the molecular topology, and reference potential energies, the models successfully learn to assign atom types and corresponding FF parameters from scratch. The resulting models and parameters are validated on a wide range of experimentally and computationally derived properties of systems including dimers, pure liquids, and molecular crystals.

Dynamics of Aqueous Electrolyte Solutions: Challenges for Simulations

This Perspective article focuses on recent simulation work on the dynamics of aqueous electrolytes. It is well-established that full-charge, nonpolarizable models for water and ions generally predict solution dynamics that are too slow in comparison to experiments. Models with reduced (scaled) charges do better for solution diffusivities and viscosities but encounter issues describing other dynamic phenomena such as nucleation rates of crystals from solution. Polarizable models show promise, especially when appropriately parametrized, but may still miss important physical effects such as charge transfer. First-principles calculations are starting to emerge for these properties that are in principle able to capture polarization, charge transfer, and chemical transformations in solution. While direct ab initio simulations are still too slow for simulations of large systems over long time scales, machine-learning models trained on appropriate first-principles data show significant promise for accurate and transferable modeling of electrolyte solution dynamics.

Accurate Hellmann-Feynman forces from density functional calculations with augmented Gaussian basis sets.

The Hellmann-Feynman (HF) theorem provides a way to compute forces directly from the electron density, enabling efficient force calculations for large systems through machine learning (ML) models for the electron density. The main issue holding back the general acceptance of the HF approach for atom-centered basis sets is the well-known Pulay force which, if naively discarded, typically constitutes an error upward of 10 eV/Å in forces. In this work, we demonstrate that if a suitably augmented Gaussian basis set is used for density functional calculations, the Pulay force can be suppressed, and HF forces can be computed as accurately as analytical forces with state-of-the-art basis sets, allowing geometry optimization and molecular dynamics to be reliably performed with HF forces. Our results pave a clear path forward for the accurate and efficient simulation of large systems using ML densities and the HF theorem.

Extending the reach of quantum computing for materials science with machine learning potentials

Solving electronic structure problems represents a promising field of applications for quantum computers. Currently, much effort is spent in devising and optimizing quantum algorithms for near-term quantum processors, with the aim of outperforming classical counterparts on selected problem instances using limited quantum resources. These methods are still expected to feature a runtime preventing quantum simulations of large scale and bulk systems. In this work, we propose a strategy to extend the scope of quantum computational methods to large scale simulations using a machine learning potential trained on quantum simulation data. The challenge of applying machine learning potentials in today’s quantum setting arises from the several sources of noise affecting the quantum computations of electronic energies and forces. We investigate the trainability of a machine learning potential selecting various sources of noise: statistical, optimization, and hardware noise. Finally, we construct the first machine learning potential from data computed on actual IBM Quantum processors for a hydrogen molecule. This already would allow us to perform arbitrarily long and stable molecular dynamics simulations, outperforming all current quantum approaches to molecular dynamics and structure optimization.

Characterization of the Interaction of Polymeric Micelles with siRNA: A Combined Experimental and Molecular Dynamics Study.

The simulation of large molecular systems remains a daunting challenge, which justifies the exploration of novel methodologies to keep computers as an ideal companion tool for everyday laboratory work. Whole micelles, bigger than 20 nm in size, formed by the self-assembly of hundreds of copolymers containing more than 50 repeating units, have until now rarely been simulated, due to a lack of computational power. Therefore, a flexible amphiphilic triblock copolymer (mPEG45-α-PLL10-PLA25) containing a total of 80 repeating units, has been emulated and synthesized to embody compactified nanoconstructs of over 900 assembled copolymers, sized between 80 and 100 nm, for siRNA complexing purposes. In this study, the tailored triblock copolymers containing a controlled number of amino groups, were used as a support model to address the binding behavior of STAT3-siRNA, in the formation of micelleplexes. Since increasingly complex drug delivery systems require an ever more optimized physicochemical characterization, a converging description has been implemented by a combination of experimentation and computational simulations. The computational data were advantageous in allowing for the assumption of an optimal N/P ratio favoring both conformational rigidifications of STAT3-siRNA with low competitive phenomena at the binding sites of the micellar carriers. These calculations were consistent with the experimental data showing that an N/P ratio of 1.5 resulted in a sufficient amount of complexed STAT3-siRNA with an electrical potential at the slipping plane of the nanopharmaceuticals, close to the charge neutralization.

Multinode Multi-GPU Two-Electron Integrals: Code Generation Using the Regent Language.

The computation of two-electron repulsion integrals (ERIs) is often the most expensive step of integral-direct self-consistent field methods. Formally it scales as O(N4), where N is the number of Gaussian basis functions used to represent the molecular wave function. In practice, this scaling can be reduced to O(N2) or less by neglecting small integrals with screening methods. The contributions of the ERIs to the Fock matrix are of Coulomb (J) and exchange (K) type and require separate algorithms to compute matrix elements efficiently. We previously implemented highly efficient GPU-accelerated J-matrix and K-matrix algorithms in the electronic structure code TeraChem. Although these implementations supported the use of multiple GPUs on a node, they did not support the use of multiple nodes. This presents a key bottleneck to cutting-edge ab initio simulations of large systems, e.g., excited state dynamics of photoactive proteins. We present our implementation of multinode multi-GPU J- and K-matrix algorithms in TeraChem using the Regent programming language. Regent directly supports distributed computation in a task-based model and can generate code for a variety of architectures, including NVIDIA GPUs. We demonstrate multinode scaling up to 45 GPUs (3 nodes) and benchmark against hand-coded TeraChem integral code. We also outline our metaprogrammed Regent implementation, which enables flexible code generation for integrals of different angular momenta.

Galvani Offset Potential and Constant-pH Simulations of Membrane Proteins.

A central problem in computational biophysics is the treatment of titratable residues in molecular dynamics simulations of large biological macromolecular systems. Conventional simulation methods ascribe a fixed ionization state to titratable residues in accordance with their pKa and the pH of the system, assuming that an effective average model will be able to capture the predominant behavior of the system. While this assumption may be justifiable in many cases, it is certainly limited, and it is important to design alternative methodologies allowing a more realistic treatment. Constant-pH simulation methods provide powerful approaches to handle titratable residues more realistically by allowing the ionization state to vary statistically during the simulation. Extending the molecular mechanical (MM) potential energy function to a family of potential functions accounting for different ionization states, constant-pH simulations are designed to sample all accessible configurations and ionization states, properly weighted according to their Boltzmann factor. Because protonation and deprotonation events correspond to a change in the total charge, difficulties arise when the long-range Coulomb interaction is treated on the basis of an idealized infinite simulation model and periodic boundary conditions with particle-mesh Ewald lattice sums. Charging free-energy calculations performed under these conditions in aqueous solution depend on the Galvani potential of the bulk water phase. This has important implications for the equilibrium and nonequilibrium constant-pH simulation methods grounded in the relative free-energy difference corresponding to the protonated and unprotonated residues. Here, the effect of the Galvani potential is clarified, and a simple practical solution is introduced to address this issue in constant-pH simulations of the acid-sensing ion channel (ASIC).

A New Transfer Free Energy Based Implicit Solvation Model for the Description of Disordered and Folded Proteins.

Most biological events occur on time scales that are difficult to access using conventional all-atom molecular dynamics simulations in explicit solvent. Implicit solvent techniques offer a promising solution to this problem, alleviating the computational cost associated with the simulation of large systems and accelerating the sampling compared to explicit solvent models. The substitution of water molecules by a mean field, however, introduces simplifications that may penalize accuracy and impede the prediction of certain physical properties. We demonstrate that existing implicit solvent models developed using a transfer free energy approach, while satisfactory at reproducing the folding behavior of globular proteins, fare less well in characterizing the conformational properties of intrinsically disordered proteins. We develop a new implicit solvent model that maximizes the degree of accuracy for both disordered and folded proteins. We show, by comparing the simulation outputs to experimental data, that in combination with the a99SB-disp force field, the implicit solvent model can describe both disordered (aβ40, PaaA2, and drkN SH3) and folded ((AAQAA)3, CLN025, Trp-cage, and GTT) peptides. Our implicit solvent model permits a computationally efficient investigation of proteins containing both ordered and disordered regions, as well as the study of the transition between ordered and disordered protein states.

Optimization of Thermal Conductance at Interfaces Using Machine Learning Algorithms.

Optimization of thermal transport across the interface of two different materials is critical to micro-/nanoscale electronic, photonic, and phononic devices. Although several examples of compositional intermixing at the interfaces having a positive effect on interfacial thermal conductance (ITC) have been reported, an optimum arrangement has not yet been determined because of the large number of potential atomic configurations and the significant computational cost of evaluation. On the other hand, computation-driven materials design efforts are rising in popularity and importance. Yet, the scalability and transferability of machine learning models remain as challenges in creating a complete pipeline for the simulation and analysis of large molecular systems. In this work we present a scalable Bayesian optimization framework, which leverages dynamic spawning of jobs through the Message Passing Interface (MPI) to run multiple parallel molecular dynamics simulations within a parent MPI job to optimize heat transfer at the silicon and aluminum (Si/Al) interface. We found a maximum of 50% increase in the ITC when introducing a two-layer intermixed region that consists of a higher percentage of Si. Because of the random nature of the intermixing, the magnitude of increase in the ITC varies. We observed that both homogeneity/heterogeneity of the intermixing and the intrinsic stochastic nature of molecular dynamics simulations account for the variance in ITC.

Spin activity correlations in driven disordered systems

We study the spatio-temporal correlations in the spin flipping activities of the disordered systems realized within the frame of the field-driven three-dimensional zero-temperature nonequilibrium random field Ising model. Our results for the spatial activity correlations reveal that the rate-dependent scaling holds in the full range of rate regimes provided that the system parameters satisfy the finite-size and rate-dependent scaling conditions. Temporal activity correlations show high sensitivity to the driving, allowing the collapsing of waiting time distributions only in the limit of slow driving rates. These findings, based on extensive numerical simulations of large systems, reliably demonstrate the observed behavior and could be relevant for studies on various, particularly experimental, nonequilibrium systems driven at finite rates.

A fast and accurate modeling approach for water and steam thermodynamics with practical applications in district heating system simulation

In U.S. district heating (DH) systems, steam is the most common heat transport medium. Industry demand for new advanced modeling capabilities of complete steam DH systems is increasing; however, the existing models for water/steam thermodynamics are too slow for large system simulations because of computationally expensive algebraic loops that require the solution to nonlinear systems of equations. For practical applications, this work presents a novel split-medium approach that implements numerically efficient liquid water models alongside various water/steam models, breaking costly algebraic loops by decoupling mass and energy balance equations. New component models for steam DH systems are also presented. We implemented the models in the equation based Modelica language and evaluated accuracy and computing speed across multiple scales: from fundamental thermodynamic properties to complete districts featuring 10 to 200 buildings. Compared to district models with the IF97 water/steam model and equipment models from the Modelica Standard Library, the new implementation improves the scaling rate for large districts from cubic to quadratic with negligible compromise to accuracy. For an annual simulation with 180 buildings, this translates to a computing time reduction from 33 to 1–1.5 h. These results are critically important for industry practitioners to simulate steam DH systems at large scales.

Path integral molecular dynamics simulations for Green's function in a system of identical bosons.

Path integral molecular dynamics (PIMD) has been successfully applied to perform simulations of large bosonic systems in a recent study [Hirshberg et al., Proc. Natl. Acad. Sci. U. S. A. 116, 21445 (2019)]. In this work, we extend PIMD techniques to study Green's function for bosonic systems. We demonstrate that the development of the original PIMD method enables us to calculate Green's function and extract momentum distribution from our simulations. We also apply our method to systems of identical interacting bosons to study Berezinskii-Kosterlitz-Thouless transition around its critical temperature.

Two-dimensional ferromagnetic systems with finite driving

We study the effect of finite driving by the time-varying external magnetic field in the two-dimensional disordered ferromagnetic systems portrayed on the case of nonequlibrium athermal random field Ising model. The response of studied systems is rate-dependent exhibiting scaling and data collapsing above the dynamical critical point in the full range of driving rates provided that the system parameters satisfy both the finite-size and here introduced rate-dependent scaling conditions. These findings are based on extensive numerical simulations of large systems reliably demonstrating the observed behavior and could be relevant for real thin systems driven at finite rates.

Self-consistent determination of long-range electrostatics in neural network potentials

Machine learning has the potential to revolutionize the field of molecular simulation through the development of efficient and accurate models of interatomic interactions. Neural networks can model interactions with the accuracy of quantum mechanics-based calculations, but with a fraction of the cost, enabling simulations of large systems over long timescales. However, implicit in the construction of neural network potentials is an assumption of locality, wherein atomic arrangements on the nanometer-scale are used to learn interatomic interactions. Because of this assumption, the resulting neural network models cannot describe long-range interactions that play critical roles in dielectric screening and chemical reactivity. Here, we address this issue by introducing the self-consistent field neural network — a general approach for learning the long-range response of molecular systems in neural network potentials that relies on a physically meaningful separation of the interatomic interactions — and demonstrate its utility by modeling liquid water with and without applied fields.

Small angle symmetry splitting of helicene-based molecular wires on pyrolytic graphite

Synthesis of inherently chiral materials and investigation of their self-assembly into molecular wires or 2D crystals belongs to hot topics in nanotechnology. Here we report an inherent splitting of orientational symmetry of 10–500 nm long molecular wires assembled from 3 nm chiral helicene-based macrocycles on atomically flat HOPG surface. The symmetry splitting of molecular wire orientation by a small angle of (6.5 ± 0.5)° is attributed to interaction of individual helicene-based macrocycle with HOPG. There is a good correlation between experimental AFM data and theoretical simulations using ReaxFF force field and Lennard-Jones (L-J) potential. In particular the new 1000 times faster computational method based on L-J potential and sequential addition of molecules is able to simulate formation of molecular wires and emergence of their symmetry splitting, which was not possible to do by conventional molecular dynamics. The method is also sensitive to chirality of the studied molecules; enantiomers and racemates placed on HOPG form different arrangements. This opens new possibilities for simulations of large molecular systems.