Simulation Of Large Systems Research Articles

Effect of three-body interactions on the zero-temperature equation of state of HCP solid 4He

Previous studies have pointed to the importance of three-body interactions in high density 4He solids. However the computational cost often makes it unfeasible to incorporate these interactions into the simulation of large systems. We report the implementation and evaluation of a computationally efficient perturbative treatment of three-body interactions in hexagonal close packed solid 4He utilizing the recently developed nonadditive three-body potential of Cencek et al. This study represents the first application of the Cencek three-body potential to condensed phase 4He systems. Ground state energies from quantum Monte Carlo simulations, with either fully incorporated or perturbatively treated three-body interactions, are calculated in systems with molar volumes ranging from 21.3 cm3/mol down to 2.5 cm3/mol. These energies are used to derive the zero-temperature equation of state for comparison against existing experimental and theoretical data. The equations of state derived from both perturbative and fully incorporated three-body interactions are found to be in very good agreement with one another, and reproduce the experimental pressure-volume data with significantly better accuracy than is obtained when only two-body interactions are considered. At molar volumes below approximately 4.0 cm3/mol, neither two-body nor three-body equations of state are able to accurately reproduce the experimental pressure-volume data, suggesting that below this molar volume four-body and higher many-body interactions are becoming important.

GPU-accelerated algorithms for many-particle continuous-time quantum walks

Many-particle continuous-time quantum walks (CTQWs) represent a resource for several tasks in quantum technology, including quantum search algorithms and universal quantum computation. In order to design and implement CTQWs in a realistic scenario, one needs effective simulation tools for Hamiltonians that take into account static noise and fluctuations in the lattice, i.e. Hamiltonians containing stochastic terms. To this aim, we suggest a parallel algorithm based on the Taylor series expansion of the evolution operator, and compare its performances with those of algorithms based on the exact diagonalization of the Hamiltonian or a 4th order Runge–Kutta integration. We prove that both Taylor-series expansion and Runge–Kutta algorithms are reliable and have a low computational cost, the Taylor-series expansion showing the additional advantage of a memory allocation not depending on the precision of calculation. Both algorithms are also highly parallelizable within the SIMT paradigm, and are thus suitable for GPGPU computing. In turn, we have benchmarked 4 NVIDIA GPUs and 3 quad-core Intel CPUs for a 2-particle system over lattices of increasing dimension, showing that the speedup provided by GPU computing, with respect to the OPENMP parallelization, lies in the range between 8x and (more than) 20x, depending on the frequency of post-processing. GPU-accelerated codes thus allow one to overcome concerns about the execution time, and make it possible simulations with many interacting particles on large lattices, with the only limit of the memory available on the device. Program summaryProgram Title: cuQuWaLicensing provisions: GNU General Public License, version 3Program Files doi:http://dx.doi.org/10.17632/vjpnjgycdj.1Programming language: CUDA CNature of problem: Evolution of many-particle continuous-time quantum-walks on a multidimensional grid in a noisy environment. The submitted code is specialized for the simulation of 2-particle quantum-walks with periodic boundary conditions.Solution method: Taylor-series expansion of the evolution operator. The density-matrix is calculated by averaging multiple independent realizations of the system.External routines: cuBLAS, cuRANDUnusual features: Simulations are run exclusively on the graphic processing unit within the CUDA environment. An undocumented misbehavior in the random-number generation routine (cuRAND package) can corrupt the simulation of large systems, though no problems are reported for small and medium-size systems. Compiling the code with the -arch=sm_30 flag for compute capability 3.5 and above fixes this issue.

Electron transport in real time from first-principles

While the vast majority of calculations reported on molecular conductance have been based on the static non-equilibrium Green's function formalism combined with density functional theory (DFT), in recent years a few time-dependent approaches to transport have started to emerge. Among these, the driven Liouville-von Neumann equation [C. G. Sánchez et al., J. Chem. Phys. 124, 214708 (2006)] is a simple and appealing route relying on a tunable rate parameter, which has been explored in the context of semi-empirical methods. In the present study, we adapt this formulation to a density functional theory framework and analyze its performance. In particular, it is implemented in an efficient all-electron DFT code with Gaussian basis functions, suitable for quantum-dynamics simulations of large molecular systems. At variance with the case of the tight-binding calculations reported in the literature, we find that now the initial perturbation to drive the system out of equilibrium plays a fundamental role in the stability of the electron dynamics. The equation of motion used in previous tight-binding implementations with massive electrodes has to be modified to produce a stable and unidirectional current during time propagation in time-dependent DFT simulations using much smaller leads. Moreover, we propose a procedure to get rid of the dependence of the current-voltage curves on the rate parameter. This method is employed to obtain the current-voltage characteristic of saturated and unsaturated hydrocarbons of different lengths, with very promising prospects.

Towards accurate quantum simulations of large systems with small computers

Numerical simulations are important for many systems. In particular, various standard computer programs have been developed for solving the quantum Schrödinger equations. However, the accuracy of these calculations is limited by computer capabilities. In this work, an iterative method is introduced to enhance the accuracy of these numerical calculations, which is otherwise prohibitive by conventional methods. The method is easily implementable and general for many systems.

Simulation of Inrush Dynamics for Unbalanced Distribution Systems Using Dynamic-Phasor Models

Microgrids are a possible solution to mitigate the extensive interruptions in service to end-use customers caused by abnormal events, such as major natural disasters. Microgrids can be formed by energizing portions of the distribution system to interconnect generators and critical loads. In this scenario, inrush dynamic currents can present a barrier for forming microgrids by preventing the necessary switching operations and/or damaging equipment. Modeling inrush dynamics, through electromagnetic transient tools, presents a challenge when the microgrid represents a large portion of a distribution feeder. This paper uses a dynamic-phasor approach to approximate the envelope of inrush dynamics; thus, allowing simulations of large systems. The two main contributions of this paper are: a) a new formulation to solve a set of dynamic-phasor differential equations representing an unbalanced distribution system; and b) modeling of transformer saturation in the dynamic-phasor frame. These two contributions allow the inrush phenomena to be simulated on full-size microgrid and/or distribution feeder models. The new formulation and model are implemented in GridLAB-D, and tested on IEEE test systems. The simulations showed good approximation as compared with a more accurate electromagnetic transient simulation in PSCAD software.

Hybrid Monte Carlo and continuum modeling of electrolytes with concentration-induced dielectric variations.

The distribution of ions near a charged surface is an important quantity in many biological and material processes, and has been therefore investigated intensively. However, few theoretical and simulation approaches have included the influence of concentration-induced variations in the local dielectric permittivity of an underlying electrolyte solution. Such local variations have long been observed and known to affect the properties of ionic solution in the bulk and around the charged surface. We propose a hybrid computational model that combines Monte Carlo simulations with continuum electrostatic modeling to investigate such properties. A key component in our hybrid model is a semianalytical formula for the ion-ion interaction energy in a dielectrically inhomogeneous environment. This formula is obtained by solving for the Green's function Poisson's equation with ionic-concentration-dependent dielectric permittivity using a harmonic interpolation method and spherical harmonic series. We also construct a self-consistent continuum model of electrostatics to describe the effect of ionic-concentration-dependent dielectric permittivity and the resulting self-energy contribution. With extensive numerical simulations, we verify the convergence of our hybrid simulation scheme, show the qualitatively different structures of ionic distribution due to the concentration-induced dielectric variations, and compare our simulation results with the self-consistent continuum model. In particular, we study the differences between weakly and strongly charged surfaces and multivalencies of counterions. Our hybrid simulations conform particularly the depletion of ionic concentrations near a charged surface and also capture the charge inversion. We discuss several issues and possible further improvement of our approach for simulations of large charged systems.

High-performance computing in systems biology: accelerating the simulation and analysis of large and complex biological systems

High-performance computing in systems biology: accelerating the simulation and analysis of large and complex biological systems

Implicit Solvent Model for Million-Atom Atomistic Simulations: Insights into the Organization of 30-nm Chromatin Fiber

Molecular dynamics (MD) simulations based on the implicit solvent generalized Born (GB) models can provide significant computational advantages over the traditional explicit solvent simulations. However, the standard GB becomes prohibitively expensive for all-atom simulations of large structures; the model scales poorly, ∼n2, with the number of solute atoms. Here we combine our recently developed optimal point charge approximation (OPCA) with the hierarchical charge partitioning (HCP) approximation to present an ∼n log n multiscale, yet fully atomistic, GB model (GB-HCPO). The HCP approximation exploits the natural organization of biomolecules (atoms, groups, chains, and complexes) to partition the structure into multiple hierarchical levels of components. OPCA approximates the charge distribution for each of these components by a small number of point charges so that the low order multipole moments of these components are optimally reproduced. The approximate charges are then used for computing electrostatic interactions with distant components, while the full set of atomic charges are used for nearby components. We show that GB-HCPO can deliver up to 2 orders of magnitude speedup compared to the standard GB, with minimal impact on its accuracy. For large structures, GB-HCPO can approach the same nominal speed, as in nanoseconds per day, as the highly optimized explicit-solvent simulation based on particle mesh Ewald (PME). The increase in the nominal simulation speed, relative to the standard GB, coupled with substantially faster sampling of conformational space, relative to the explicit solvent, makes GB-HCPO a suitable candidate for MD simulation of large atomistic systems in implicit solvent. As a practical demonstration, we use GB-HCPO simulation to refine a ∼1.16 million atom structure of 30 nm chromatin fiber (40 nucleosomes). The refined structure suggests important details about spatial organization of the linker DNA and the histone tails in the fiber: (1) the linker DNA fills the core region, allowing the H3 histone tails to interact with the linker DNA, which is consistent with experiment; (2) H3 and H4 tails are found mostly in the core of the structure, closer to the helical axis of the fiber, while H2A and H2B are mostly solvent exposed. Potential functional consequences of these findings are discussed. GB-HCPO is implemented in the open source MD software NAB in Amber 2016.

DNS–Assisted Modeling of Bubbly Flows in Vertical Channels

The transient motion of bubbly flows in vertical channels is studied, using direct numerical simulation (DNS) in which every continuum length and time scale is resolved. A simulation of a large number of bubbles of different sizes at a friction Reynolds number of 500 shows that small bubbles quickly migrate to the wall, but the bulk flow takes much longer to adjust to the new bubble distribution. Simulations of much smaller laminar systems with several spherical bubbles have been used to examine the full transient motion; those show a nonmonotonic evolution where all the bubbles first move toward the walls, and the liquid then slowly slows down, eventually allowing some bubbles to return to the center of the channel. Unlike the statistically steady state, where the flow structure is relatively simple and in some cases depends only on the sign of the bubble lift coefficient, the transient evolution is more sensitive to the governing parameters. Early efforts to use DNS results to provide values for the unresolved closure terms in a simple average model for the flow found by statistical learning from the data using neural networks are discussed. The prospect for using the results from simulations of large systems with bubbles of different sizes in turbulent flows for large eddy—like simulations are explored, including the simplification of the interface structure by filtering. Finally, preliminary results for flows undergoing topology changes are shown.

Problems and Development of Methods of Dynamic Simulation of Cryogenic Systems

The problems of designing and employing cryogenic systems operating in transient modes over a considerable operating time as well as various methods of simulating transient operation of these systems are discussed. The advantages of application of dynamic simulation of cryogenic systems over classical methods of calculation of steady-state mode are shown. The available information on dynamic simulation of large helium cryogenic systems is reviewed. To validate the dynamic simulation method, the transient mode of operation of a large helium liquefier is analyzed applying universal software tool for dynamic simulation modeling of chemical systems, the results of which are compared with the experimental data.

Graphics Processing Unit Acceleration and Parallelization of GENESIS for Large-Scale Molecular Dynamics Simulations.

The graphics processing unit (GPU) has become a popular computational platform for molecular dynamics (MD) simulations of biomolecules. A significant speedup in the simulations of small- or medium-size systems using only a few computer nodes with a single or multiple GPUs has been reported. Because of GPU memory limitation and slow communication between GPUs on different computer nodes, it is not straightforward to accelerate MD simulations of large biological systems that contain a few million or more atoms on massively parallel supercomputers with GPUs. In this study, we develop a new scheme in our MD software, GENESIS, to reduce the total computational time on such computers. Computationally intensive real-space nonbonded interactions are computed mainly on GPUs in the scheme, while less intensive bonded interactions and communication-intensive reciprocal-space interactions are performed on CPUs. On the basis of the midpoint cell method as a domain decomposition scheme, we invent the single particle interaction list for reducing the GPU memory usage. Since total computational time is limited by the reciprocal-space computation, we utilize the RESPA multiple time-step integration and reduce the CPU resting time by assigning a subset of nonbonded interactions on CPUs as well as on GPUs when the reciprocal-space computation is skipped. We validated our GPU implementations in GENESIS on BPTI and a membrane protein, porin, by MD simulations and an alanine-tripeptide by REMD simulations. Benchmark calculations on TSUBAME supercomputer showed that an MD simulation of a million atoms system was scalable up to 256 computer nodes with GPUs.

A fast integral equation method for solid particles in viscous flow using quadrature by expansion

Boundary integral methods are advantageous when simulating viscous flow around rigid particles, due to the reduction in number of unknowns and straightforward handling of the geometry. In this work we present a fast and accurate framework for simulating spheroids in periodic Stokes flow, which is based on the completed double layer boundary integral formulation. The framework implements a new method known as quadrature by expansion (QBX), which uses surrogate local expansions of the layer potential to evaluate it to very high accuracy both on and off the particle surfaces. This quadrature method is accelerated through a newly developed precomputation scheme. The long range interactions are computed using the spectral Ewald (SE) fast summation method, which after integration with QBX allows the resulting system to be solved in Mlog⁡M time, where M is the number of particles. This framework is suitable for simulations of large particle systems, and can be used for studying e.g. porous media models.

Phase transitions and ordering structures of a model of a chiral helimagnet in three dimensions

Phase transitions in a classical Heisenberg spin model of a chiral helimagnet with the Dzyaloshinskii--Moriya (DM) interaction in three dimensions are numerically studied. By using the event-chain Monte Carlo algorithm recently developed for particle and continuous spin systems, we perform equilibrium Monte Carlo simulations for large systems up to about $10^6$ spins. Without magnetic fields, the system undergoes a continuous phase transition with critical exponents of the three-dimensional \textit{XY} model, and a uniaxial periodic helical structure emerges in the low temperature region. In the presence of a magnetic field perpendicular to the axis of the helical structure, it is found that there exists a critical point on the temperature and magnetic-field phase diagram and that above the critical point the system exhibits a phase transition with strong divergence of the specific heat and the uniform magnetic susceptibility.

Experimental investigation and modeling of a hermetic scroll expander

As part of the Enerco_LT (ENErgy ReCOvery from Low Temperature heat sources) project, hermetic scroll expander is candidate for small-scale ORC prototypes. The aim of this work is (i) to investigate on a test bench the performances of a small-size scroll expander which is a converted compressor and (ii) to develop, based on the experimental results, a semi-empirical steady-state model which can be easily integrated as a sub-model in a larger system simulation model.The expander was tested by means of a superheated vapour cycle (Brayton) fed with R-245fa and providing defined ranges of supply and exhaust operating conditions. The assessment of the internal unknown variables has permitted the derivation of the model correlations. The supply pressure drop and the internal isentropic efficiency are expressed as polynomial laws of the main working conditions. The model was validated by comparing predictions to the experimental results which has shown a fairly good agreement and accuracy.The experimental investigation has revealed promising performances of the hermetic scroll expander given at high supply conditions. The net efficiency of the combined expander-generator reaches 70% from supply pressures of 20bar. Furthermore, the numerical investigation using the developed model has pointed out the performance deterioration when the pressure ratio deviates from the expander built-in volume ratio and indicates consequently how optimizing such machines to ensure the highest performances.

High-Efficiency Microiterative Optimization in QM/MM Simulations of Large Flexible Systems.

We present here a double-optimizations-of-buffer-region (DOBR) microiterative scheme for high-efficiency energy minimizations of large, flexible systems in combined quantum-mechanical/molecular-mechanical (QM/MM) calculations. In the DOBR scheme, an entire system is divided into three regions: the core, buffer, and outer regions. The core region includes QM atoms and the MM atoms within a cutoff distance R1 to the QM atoms (denoted by MM1 atoms), and the buffer region consists of MM atoms within another cutoff distance R2 to MM1 atoms. Each DOBR microcycle involves two steps: First, QM atoms are assigned electrostatic-potential (ESP) charges, and the buffer and outer regions are optimized at the MM level with the core region kept frozen. Second, the core and buffer regions are optimized at the QM/MM level using the electrostatic embedding with the outer region kept frozen. The two steps are repeated until two optimizations converge at one structure. The DOBR scheme was tested in the optimizations of nucleobases solvated in water spheres of 30 Å radius, where the initial geometries were extracted from the trajectories of classical molecular dynamics simulations, and the cutoff distances R1 and R2 were set to 5.0 and 4.0 Å, respectively. For comparisons, the optimizations were also carried out by a "standard" scheme without microiteration and by the two-region microiterative (TRM) method. We found that the averaged number of QM calculations for the DOBR scheme is only ∼1% of that of the standard scheme and ∼6% of the TRM approach. The promising results indicate that the DOBR scheme could significantly increase the efficiency of geometry optimizations for large, flexible systems in QM/MM calculations.

ClustENM: ENM-Based Sampling of Essential Conformational Space at Full Atomic Resolution.

Accurate sampling of conformational space and, in particular, the transitions between functional substates has been a challenge in molecular dynamic (MD) simulations of large biomolecular systems. We developed an Elastic Network Model (ENM)-based computational method, ClustENM, for sampling large conformational changes of biomolecules with various sizes and oligomerization states. ClustENM is an iterative method that combines ENM with energy minimization and clustering steps. It is an unbiased technique, which requires only an initial structure as input, and no information about the target conformation. To test the performance of ClustENM, we applied it to six biomolecular systems: adenylate kinase (AK), calmodulin, p38 MAP kinase, HIV-1 reverse transcriptase (RT), triosephosphate isomerase (TIM), and the 70S ribosomal complex. The generated ensembles of conformers determined at atomic resolution show good agreement with experimental data (979 structures resolved by X-ray and/or NMR) and encompass the subspaces covered in independent MD simulations for TIM, p38, and RT. ClustENM emerges as a computationally efficient tool for characterizing the conformational space of large systems at atomic detail, in addition to generating a representative ensemble of conformers that can be advantageously used in simulating substrate/ligand-binding events.

Many-body localization and transition by density matrix renormalization group and exact diagonalization studies

A many-body localized (MBL) state is a new state of matter emerging in a disordered interacting system at high energy densities through a disorder driven dynamic phase transition. The nature of the phase transition and the evolution of the MBL phase near the transition are the focus of intense theoretical studies with open issues in the field. We develop an entanglement density matrix renormalization group (En-DMRG) algorithm to accurately target the entanglement patterns of highly excited states for MBL systems. By studying the one dimensional Heisenberg spin chain in a random field, we demonstrate the high accuracy of the method in obtaining statistical results of quantum states in the MBL phase. Based on large system simulations by En-DMRG for excited states, we demonstrate some interesting features in the entanglement entropy distribution function, which is characterized by two peaks; one at zero and another one at the quantized entropy $S=ln2$ with an exponential decay tail on the $S>ln2$ side. Combining En-DMRG with exact diagonalization simulations, we demonstrate that the transition from the MBL phase to the delocalized ergodic phase is driven by rare events where the locally entangled spin pairs develop long-range power-law correlations. The corresponding phase diagram contains an intermediate regime, which has power-law spin-z correlations resulting from contributions of the rare events. We discuss the physical picture for the numerical observations in the intermediate regime, where various distribution functions are distinctly different from results deep in the ergodic and MBL phases for finite-size systems. Our results may provide new insights for understanding the phase transition in such systems.

A QSAR model of benzoxazole derivatives as potential inhibitors for inosine 5`-monophosphate dehydrogenase from Cryptosporidium parvum.

Cryptosporidium parvum is the common enteric protozoan pathogen causing cryptosporidiosis in human. Available drugs to treat cryptosporidiosis are ineffective and there is yet no vaccine against C. parvum. Therefore, it is of interest to design an improved yet effective drug against C. parvum. Here, we docked benzoxazole derivatives (collected from literature) with inosine 5`- monophosphate dehydrogenase (IMPDH) from Cryptosporidium parvum using the program AutoDock 4.2. The docked protein - inhibitor complex structure was optimized using molecular dynamics simulation for 5 ps with the CHARMM-22 force field using NAMD (NAnoscale Molecular Dynamics program) incorporated in visual molecular dynamics (VMD 1.9.2) and then evaluating the stability of complex structure by calculating RMSD values. NAMD is a parallel, object-oriented molecular dynamics code designed for high-performance simulation of large biomolecular systems. A quantitative structure activity relationship (QSAR) model was built using energy-based descriptors as independent variable and pIC50 value as dependent variable of fifteen known benzoxazole derivatives with C. parvum IMPDH protein, yielding correlation coefficient r2 of 0.7948. The predictive performance of QSAR model was assessed using different cross-validation procedures. Our results suggest that a ligand-receptor binding interaction for inosine 5`-monophosphate dehydrogenase using a QSAR model is promising approach to design more potent inosine 5`-monophosphate dehydrogenase inhibitors prior to their synthesis.

Combining Linear-Scaling DFT with Subsystem DFT in Born-Oppenheimer and Ehrenfest Molecular Dynamics Simulations: From Molecules to a Virus in Solution.

In this work, methods for the efficient simulation of large systems embedded in a molecular environment are presented. These methods combine linear-scaling (LS) Kohn-Sham (KS) density functional theory (DFT) with subsystem (SS) DFT. LS DFT is efficient for large subsystems, while SS DFT is linear scaling with a smaller prefactor for large sets of small molecules. The combination of SS and LS, which is an embedding approach, can result in a 10-fold speedup over a pure LS simulation for large systems in aqueous solution. In addition to a ground-state Born-Oppenheimer SS+LS implementation, a time-dependent density functional theory-based Ehrenfest molecular dynamics (EMD) using density matrix propagation is presented that allows for performing nonadiabatic dynamics. Density matrix-based EMD in the SS framework is naturally linear scaling and appears suitable to study the electronic dynamics of molecules in solution. In the LS framework, linear scaling results as long as the density matrix remains sparse during time propagation. However, we generally find a less than exponential decay of the density matrix after a sufficiently long EMD run, preventing LS EMD simulations with arbitrary accuracy. The methods are tested on various systems, including spectroscopy on dyes, the electronic structure of TiO2 nanoparticles, electronic transport in carbon nanotubes, and the satellite tobacco mosaic virus in explicit solution.

Effect of van der Waals force cut-off distance on adhesive collision parameters in DEM simulation

Rapid advancement in computer technology makes it possible to perform simulations of large particulate systems using Discrete Element Method (DEM). However, it is still a challenge to adjust the DEM simulation parameters, especially the interparticle forces, in order to obtain rational results. Many of these forces diverge as the distance between particle surfaces approaches zero. Therefore, a cut-off distance needs to be considered in order to avoid such a problem. This work presents the results of the computational analysis of the influence of the cut-off distance, which is required to avoid the divergence of the van der Waals force, on the collision of a single particle against a flat surface. Hence, the effect of the cut-off distance on the coefficient of restitution, collision duration and maximum overlap has been studied. In addition, we have theoretically derived an expression for the minimum velocity under which the particle remains adhered to the surface (critical velocity) as a function of the cut-off distance. The simulation predictions of the critical velocity are within the range of experimental data published in the literature. We demonstrate that the cut-off distance has a profound influence on particle rebound and therefore, a careful selection of this parameter should take place when simulating bulk particle behaviour. Given that the hydrophobic force is usually simulated using the same expression as the van der Waals force, the results presented here could also be considered in the context of the simulation of the hydrophobic force.