Efficient electronic structure theory via hierarchical scale-adaptive coupled-cluster formalism: I. Theory and computational complexity analysis

ABSTRACT

A novel reduced-scaling, general-order coupled-cluster approach is formulated by exploiting hierarchical representations of many-body tensors, combined with the recently suggested formalism of scale-adaptive tensor algebra. Inspired by the hierarchical techniques from the renormalisation group approach, H/H²-matrix algebra and fast multipole method, the computational scaling reduction in our formalism is achieved via coarsening of quantum many-body interactions at larger interaction scales, thus imposing a hierarchical structure on many-body tensors of coupled-cluster theory. In our approach, the interaction scale can be defined on any appropriate Euclidean domain (spatial domain, momentum-space domain, energy domain, etc.). We show that the hierarchically resolved many-body tensors can reduce the storage requirements to O(N), where N is the number of simulated quantum particles. Subsequently, we prove that any connected many-body diagram consisting of a finite number of arbitrary-order tensors, e.g. an arbitrary coupled-cluster diagram, can be evaluated in O(NlogN) floating-point operations. On top of that, we suggest an additional approximation to further reduce the computational complexity of higher order coupled-cluster equations, i.e. equations involving higher than double excitations, which otherwise would introduce a large prefactor into formal O(NlogN) scaling.     

A subgrid model for clustering of high-inertia particles in large-eddy simulations of turbulence

Clustering (or preferential concentration) of inertial particles suspended in a homogeneous, isotropic turbulent flow is strongly influenced by the smallest scales of the turbulence. In particle-laden large-eddy simulations (LES) of turbulence, these small scales are not captured by the grid and hence their effect on particle motion needs to be modelled. In this paper, we use a subgrid model based on kinematic simulations of turbulence (Kinematic Simulation based SubGrid Model or KSSGM), for the first time in the context of predicting the clustering and the relative velocity statistics of inertial particles. This initial study focuses on the special case of inertial particles in the absence of gravitational settling. We show that the KSSGM gives excellent predictions for clustering in a priori tests for inertial particles with St ≥ 2.0, where St is the Stokes number, defined as the ratio of the particle response time to the Kolmogorov time-scale. To the best of our knowledge, the KSSGM represents the first model that has been shown to capture the effect of the subgrid scales on inertial particle clustering for St ≥ 2.0. We also show that the mean inward radial relative velocity between inertial particles (?w_r?^(?), which enters into the formula for the collision kernel) is accurately predicted by the KSSGM for all St. We explain why the model captures clustering at higher St?but not for lower St?, and provide new insights into the key statistical parameters of turbulence that a subgrid model would have to describe, in order to accurately predict clustering of low-St?particles in an LES.     

Towardsab Initio simulations of concentrated suspensions

Determination of many-body interactions between particles of arbitrary shape in a viscous fluid is a key problem in the simulation of concentrated suspensions. Three-dimensional flows involving such complex fluid-solid boundaries are beyond the scope of spatial methods, even on supercomputers. Boundary integral methods convert the three-dimensional PDE to a two-dimensional integral equation. Unfortunately, conventional boundary methods yield Fredholm integral equations of the first kind, and dense linear systems which are too large for accurate solution. We have pursued a different boundary integral formulation, which yields Fredholm integral equations of the second kind; these arc amenable to iterative solution. The velocity representation involves a compact operator, so a discrete spectrum results. Wielandt deflations give dramatic reductions in the spectral radius and accurate solutions are obtained after only a few iterations (typically less than 10). An analytic construction of the spectrum for sphere sphere interactions confirms these numerical results. The mathematics is similar to that encountered in the mixing ofd-atomic orbitals to form bonding/antibonding molecular orbitals in transition metals. The memory-saving version of our code can be implemented directly on a dedicated MicroVAX to solve problems involving clusters of less than a dozen particles. For a fixed number of processors, the algorithm grows essentially asN ², whereN is the system size, so computational times are readily estimated on more powerful super-minicomputers and supercomputers using standard dot-product benchmarks. The algorithm is especially ideal for gigaflop and teraflop parallel array processors under construction in a number of computer companies; an analysis of the spectrum reveals that asynchronous iterative methods will converge, leading the way to a rigorous formulation of screening concepts for suspended particles of arbitrary shape.     

Species segregation in one-dimensional granular-system simulations

We present one-dimensional molecular dynamics simulations of a two-species, initially uniform, freely evolving granular system. Colliding particles swap their relative position with a 50% probability allowing for the initial spatial ordering of the particles to evolve in time and frictional forces to operate. Unlike one-dimensional systems of identical particles, two-species one-dimensional systems of quasi-elastic particles are ergodic and the particles' velocity distributions tend to evolve towards Maxwell-Boltzmann distributions. Under such conditions, standard fluid equations with merely an additional sink term in the energy equation, reflecting the non-elasticity of the interparticle collisions, provide an excellent means to investigate the system's evolution. According to the predictions of fluid theory we find that the clustering instability is dominated by a non-propagating mode at a wavelength of the order 10πL/Nɛ , where N is the total number of particles, L the spatial extent of the system and ɛ the inelasticity coefficient. The typical fluid velocities at the time of inelastic collapse are seen to be supersonic, unless Nɛ ≲ 10π . Species segregation, driven by the frictional force occurs as a result of the strong temperature gradients within clusters which pushes the light particles towards the clusters' edges and the heavy particles towards the center. Segregation within clusters is complete at the time of inelastic collapse.     

Modeling of clusters by a molecular dynamics model using a fast tree method

The dynamics of clusters irradiated by a high-intensity ultrashort pulse laser has been studied using a fully relativistic three-dimensional Molecular Dynamics Model. A fast three-dimensional tree algorithm for computing the electrostatic force has been developed and compared with the conventional particle-particle method. The particle-particle method requires computation time, which scales as O(N_p ²), and it is faster for small number of particles N_p <10³. In the opposite case of relatively large ensemble of particles N_p >10³, the preferred method is the tree algorithm whose computation time scales as O(N_p log N_p). The tree algorithm has been benchmarked against the particle-particle method for clusters composed of xenon and deuterium atoms and its accuracy and computation time have been analyzed. The optimum free parameter of the tree method has been determined to be θ≈0.5. We addressed the effects of boundary conditions by studying the contribution of adjacent clusters to the total electromagnetic force exerted on individual particles. We found that the adjacent clusters play a minor role in the overall cluster dynamics.     

Critical anomaly and finite size scaling of the self-diffusion coefficient for Lennardben Jones fluids by non-equilibrium molecular dynamic simulation

We use non-equilibrium molecular dynamics simulations to calculate the self-diffusion coefficient, D, of a Lennard-Jones fluid over a wide density and temperature range. The change in self-diffusion coefficient with temperature decreases by increasing density. For density ρ^* = ρσ³ = 0.84 we observe a peak at the value of the self-diffusion coefficient and the critical temperature T^* = kT/ε = 1.25. The value of the self-diffusion coefficient strongly depends on system size. The data of the self-diffusion coefficient are fitted to a simple analytic relation based on hydrodynamic arguments. This correction scales as N^-α, where α is an adjustable parameter and N is the number of particles. It is observed that the values of α < 1 provide quite a good correction to the simulation data. The system size dependence is very strong for lower densities, but it is not as strong for higher densities. The self-diffusion coefficient calculated with non-equilibrium molecular dynamic simulations at different temperatures and densities is in good agreement with other calculations from the literature.     

Euler's fluid equations: Optimal control vs optimization

An optimization method used in image-processing (metamorphosis) is found to imply Euler's equations for incompressible flow of an inviscid fluid, without requiring that the Lagrangian particle labels exactly follow the flow lines of the Eulerian velocity vector field. Thus, an optimal control problem and an optimization problem for incompressible ideal fluid flow both yield the same Euler fluid equations, although their Lagrangian parcel dynamics are different. This is a result of the gauge freedom in the definition of the fluid pressure for an incompressible flow, in combination with the symmetry of fluid dynamics under relabeling of their Lagrangian coordinates. Similar ideas are also illustrated for SO(N) rigid body motion.     

Links between microscopic and macroscopic fluid mechanics

The microscopic and macroscopic versions of fluid mechanics differ qualitatively. Microscopic particles obey time-reversible ordinary differential equations. The resulting particle trajectories {q(t)} may be time-averaged or ensemble-averaged so as to generate field quantities corresponding to macroscopic variables. On the other hand, the macroscopic continuum fields described by fluid mechanics follow irreversible partial differential equations. Smooth particle methods bridge the gap separating these two views of fluids by solving the macroscopic field equations with particle dynamics that resemble molecular dynamics. Recently, nonlinear dynamics have provided some useful tools for understanding the relationship between the microscopic and macroscopic points of view. Chaos and fractals play key roles in this new understanding. Non-equilibrium phase-space averages look very different from their equilibrium counterparts. Away from equilibrium the smooth phase-space distributions are replaced by fractional-dimensional singular distributions that exhibit time irreversibility.     

Identification of the Dilute Regime in Particle Sedimentation

We investigate the dynamics of rigid, spherical particles of radius R sinking in a viscous fluid. Both the inertia of the particles and the fluid are neglected. We are interested in a large number N of particles with average distance dR. We investigate in which regime (in terms of N and R/d) the particles do not significantly interact and approximately sink like single particles. We rigorously establish the lower bound for the critical number N_crit of particles. This lower bound agrees with the heuristically expected N_crit in terms of its scaling in R/d. The main difficulty lies in showing that the particles cannot get significantly closer over a relevant time scale. We use the method of reflection for the Stokes operator to bound the strength of the particle interaction.     

Crossover from normal diffusion to single-file diffusion of particles in a one-dimensional channel: LJ particles in zeolite zsm-22

The crossover from normal Fickian diffusion, where mean-squared displacement goes as t to single-file diffusion, where <z²> t?^0.5 is studied as function of particle size confined in zeolite zsm-22 using molecular dynamics simulations. The simulation results indicate that the crossover is smooth as the particle size increases and the diffusion reaches single file through a series of subdiffusion processes. A small number of mutual passage of particles can destroy the correlated single-file diffusion. The density dependence on single-file mobility obtained from molecular dynamics simulations are in very good agreement with theoretical predictions.     

Fokker–Planck Equation,Molecular Friction,and Molecular Dynamics for Brownian Particle Transport near External Solid Surfaces

The Fokker–Planck (FP) equation describing the dynamics of a single Brownian particle near a fixed external surface is derived using the multiple-time-scales perturbation method, previously used by Cukier and Deutch and Nienhuis in the absence of any external surfaces, and Piasecki et al. for two Brownian spheres in a hard fluid. The FP equation includes an explicit expression for the (time-independent) particle friction tensor in terms of the force autocorrelation function and equilibrium average force on the particle by the surrounding fluid and in the presence of a fixed external surface, such as an adsorbate. The scaling and perturbation analysis given here also shows that the force autocorrelation function must decay rapidly on the zeroth-order time scale ₀, which physically requires N _Kn1, where N _Kn is the Knudsen number (ratio of the length scale for fluid intermolecular interactions to the Brownian particle length scale). This restricts the theory given here to liquid systems where N _Kn1. For a specified particle configuration with respect to the external surface, equilibrium canonical molecular dynamics (MD) calculations are conducted, as shown here, in order to obtain numerical values of the friction tensor from the force autocorrelation expression. Molecular dynamics computations of the friction tensor for a single spherical particle in the absence of a fixed external surface are shown to recover Stokes' law for various types of fluid molecule–particle interaction potentials. Analytical studies of the static force correlation function also demonstrate the remarkable principle of force-time parity whereby the particle friction coefficient is nearly independent of the fluid molecule–particle interaction potential. Molecular dynamics computations of the friction tensor for a single spherical particle near a fixed external spherical surface (adsorbate) demonstrate a breakdown in continuum hydrodynamic results at close particle–surface separation distances on the order of several molecular diameters.     

An O(N logN) alternating-direction finite difference method for two-dimensional fractional diffusion equations

Fractional diffusion equations model phenomena exhibiting anomalous diffusion that cannot be modeled accurately by the second-order diffusion equations. Because of the nonlocal property of fractional differential operators, the numerical methods for fractional diffusion equations often generate dense or even full coefficient matrices. Consequently, the numerical solution of these methods often require computational work of O(N³) per time step and memory of O(N²) for where N is the number of grid points.     

Dynamic adaptive chemistry for turbulent flame simulations

The use of large chemical mechanisms in flame simulations is computationally expensive due to the large number of chemical species and the wide range of chemical time scales involved. This study investigates the use of dynamic adaptive chemistry (DAC) for efficient chemistry calculations in turbulent flame simulations. DAC is achieved through the directed relation graph (DRG) method, which is invoked for each computational fluid dynamics cell/particle to obtain a small skeletal mechanism that is valid for the local thermochemical condition. Consequently, during reaction fractional steps, one needs to solve a smaller set of ordinary differential equations governing chemical kinetics. Test calculations are performed in a partially-stirred reactor (PaSR) involving both methane/air premixed and non-premixed combustion with chemistry described by the 53-species GRI-Mech 3.0 mechanism and the 129-species USC-Mech II mechanism augmented with recently updated NO _x pathways, respectively. Results show that, in the DAC approach, the DRG reduction threshold effectively controls the incurred errors in the predicted temperature and species concentrations. The computational saving achieved by DAC increases with the size of chemical kinetic mechanisms. For the PaSR simulations, DAC achieves a speedup factor of up to three for GRI-Mech 3.0 and up to six for USC-Mech II in simulation time, while at the same time maintaining good accuracy in temperature and species concentration predictions.     

A guide to RBF-generated finite differences for nonlinear transport: Shallow water simulations on a sphere

The current paper establishes the computational efficiency and accuracy of the RBF-FD method for large-scale geoscience modeling with comparisons to state-of-the-art methods as high-order discontinuous Galerkin and spherical harmonics, the latter using expansions with close to 300,000 bases. The test cases are demanding fluid flow problems on the sphere that exhibit numerical challenges, such as Gibbs phenomena, sharp gradients, and complex vortical dynamics with rapid energy transfer from large to small scales over short time periods. The computations were possible as well as very competitive due to the implementation of hyperviscosity on large RBF stencil sizes (corresponding roughly to 6th to 9th order methods) with up to O(10⁵) nodes on the sphere. The RBF-FD method scaled as O(N) per time step, where N is the total number of nodes on the sphere. In Appendix A, guidelines are given on how to chose parameters when using RBF-FD to solve hyperbolic PDEs.     

The Mean-Field Limit for Solid Particles in a Navier-Stokes Flow

We propose a mathematical derivation of Brinkman’s force for a cloud of particles immersed in an incompressible viscous fluid. Specifically, we consider the Stokes or steady Navier-Stokes equations in a bounded domain Ω⊂ℝ³ for the velocity field u of an incompressible fluid with kinematic viscosity ν and density 1. Brinkman’s force consists of a source term 6π ν j where j is the current density of the particles, and of a friction term 6π ν ρ u where ρ is the number density of particles. These additional terms in the motion equation for the fluid are obtained from the Stokes or steady Navier-Stokes equations set in Ω minus the disjoint union of N balls of radius ε=1/N in the large N limit with no-slip boundary condition. The number density ρ and current density j are obtained from the limiting phase space empirical measure , where x _k is the center of the k-th ball and v _k its instantaneous velocity. This can be seen as a generalization of Allaire’s result in [Arch. Ration. Mech. Anal. 113:209–259, [1991]] who considered the case of periodically distributed x _k s with v _k =0, and our proof is based on slightly simpler though similar homogenization arguments. Similar equations are used for describing the fluid phase in various models for sprays.     

A direct O(N log N) finite difference method for fractional diffusion equations

Fractional diffusion equations model phenomena exhibiting anomalous diffusion that can not be modeled accurately by the second-order diffusion equations. Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require storage of O(N²) and computational cost of O(N³) where N is the number of grid points.     

Locating the minimum: Approach to equilibrium in a disordered, symmetric zero range process

We consider the dynamics of the disordered, one-dimensional, symmetric zero range process in which a particle from an occupied site k hops to its nearest neighbor with a quenched rate w(k). These rates are chosen randomly from the probability distribution f(w) ∼ (w − c) ⁿ , where c is the lower cutoff. For n>0, this model is known to exhibit a phase transition in the steady state from a low density phase with a finite number of particles at each site to a high density aggregate phase in which the site with the lowest hopping rate supports an infinite number of particles. In the latter case, it is interesting to ask how the system locates the site with globally minimum rate. We use an argument based on the local equilibrium, supported by Monte Carlo simulations, to describe the approach to the steady state. We find that at large enough time, regions with a smooth density profile are described by a diffusion equation with site-dependent rates, while the isolated points where the mass distribution is singular act as the boundaries of these regions. Our argument implies that the relaxation time scales with the system size L as L ^z with z = 2 + 1/(n + 1) for n>1 and suggests a different behavior for n<1.     

Exponential Convergence to the Maxwell Distribution for Some Class of Boltzmann Equations

We consider a class of nonlinear Boltzmann equations describing return to thermal equilibrium in a gas of colliding particles suspended in a thermal medium. We study solutions in the space where is the one-particle phase space and is the Liouville measure on Γ⁽¹⁾. Special solutions of these equations, called “Maxwellians,” are spatially homogenous static Maxwell velocity distributions at the temperature of the medium. We prove that, for dilute gases, the solutions corresponding to smooth initial conditions in a weighted L ¹-space converge to a Maxwellian in , exponentially fast in time.     

Crossover from Quantum to Boltzmann Statistics for Free Particles in a Single Harmonic Trap

With invoking analytical formulae in number theory and numerical calculations, we calculate the number of microstates in microcanonical ensemble for free particles in a single harmonic trap which in whole space defines a thermodynamic system but not a spatially homogeneous one. Once the number of excitation quanta m is larger than the square of the particle number N ² as m⪰O(N ²) when N≫1, the number of microcanonical microstates for an ideal, harmonically trapped Bose or Fermi gas gradually converge to the Boltzmann microcanonical microstates for the classical particles with a proper consideration of the indistinguishability.