Search

We present for the first time an efficient iterative method to directly solve the four-component Dirac-Kohn-Sham (DKS) density functional theory. Due to the existence of the negative energy continuum in the DKS operator, the existing iterative techniques for solving the Kohn-Sham systems cannot be efficiently applied to solve the DKS systems. The key component of our method is a novel filtering step (F) which acts as a preconditioner in the framework of the locally optimal block preconditioned conjugate gradient (LOBPCG) method. The resulting method, dubbed the LOBPCG-F method, is able to compute the desired eigenvalues and eigenvectors in the positive energy band without computing any state in the negative energy band. The LOBPCG-F method introduces mild extra cost compared to the standard LOBPCG method and can be easily implemented. We demonstrate our method in the pseudopotential framework with a planewave basis set which naturally satisfies the kinetic balance prescription. Numerical results for Pt$_{2}$, Au$_{2}$, TlF, and Bi$_{2}$Se$_{3}$ indicate that the LOBPCG-F method is a robust and efficient method for investigating the relativistic effect in systems containing heavy elements., 31 pages, 5 figures

We develop a technique for generating a set of optimized local basis functions to solve models in the Kohn-Sham density functional theory for both insulating and metallic systems. The optimized local basis functions are obtained by solving a minimization problem in an admissible set determined by a large number of primitive basis functions. Using the optimized local basis set, the electron energy and the atomic force can be calculated accurately with a small number of basis functions. The Pulay force is systematically controlled and is not required to be calculated, which makes the optimized local basis set an ideal tool for ab initio molecular dynamics and structure optimization. We also propose a preconditioned Newton-GMRES method to obtain the optimized local basis functions in practice. The optimized local basis set is able to achieve high accuracy with a small number of basis functions per atom when applied to a one dimensional model problem.

We present a multiscale model for numerical simulations of dynamics of crystalline solids. The method combines the continuum nonlinear elasto-dynamics model, which models the stress waves and physical loading conditions, and molecular dynamics model, which provides the nonlinear constitutive relation and resolves the atomic structures near local defects. The coupling of the two models is achieved based on a general framework for multiscale modeling - the heterogeneous multiscale method (HMM). We derive an explicit coupling condition at the atomistic/continuum interface. Application to the dynamics of brittle cracks under various loading conditions is presented as test examples.

A numerical approach based on the string method is developed to study nucleation of ordered phases in first-order phase transitions. Among other things, this method allows an efficient computation of the minimum energy path (MEP) during the nucleation process. The MEP provides information about the size, shape and free energy barrier of the critical nucleus. To improve the efficiency of the string method, a special initialization process is proposed. Constraints from physical models are treated using two methods, a generalized coordinates method and a projection method. Strategies for choosing the computational domain and defining the nucleus boundary are also introduced. The validity of our approach is illustrated by two nontrivial examples from soft condensed matter physics, namely the nematic-isotropic transition of liquid crystals and the ordered-to-ordered phase transition of diblock copolymers.

We present a new general framework for designing multiscale methods. Compared with previous work such as Brandt's systematic up-scaling, the heterogeneous multiscale method (HMM) and the ''equation-free'' approach, this new framework has the distinct feature that it does not require reinitializing the microscale model at each macro time step or each macro iteration step. In the new strategy, the macro- and micro-models evolve simultaneously using different time steps (and therefore different clocks), and they exchange data at every step. The micro-model uses its own appropriate time step. The macro-model runs at a slower pace than required by accuracy and stability considerations for the macroscale dynamics, in order for the micro-model to relax. Examples are discussed and application to modeling complex fluids is presented.

In this paper we develop a domain decomposition method (DDM), based on the discontinuous Galerkin (DG) and the local discontinuous Galerkin (LDG) methods, for solving multiscale problems involving macro sub-domains, where a macro model is valid, and micro sub-domains, where the macro model is not valid and a more costly micro model must be used. We take two examples, one from compressible gas dynamics where the micro sub-domains are around shocks, contacts and corners of rarefaction fans, and another one from semiconductor device simulations where the micro sub-domains are around the jumps in the doping profile. The macro model is taken as the Euler equations for the gas dynamics problem and as a hydrodynamic model and a high field model for the semiconductor device problem. The micro model for both problems is taken as a kinetic equation. We pay special attention to the effective coupling between the macro sub-domains and the micro sub-domains, in which we utilize the advantage of the discontinuous Galerkin method in its compactness of the computational stencil. Numerical results demonstrate the effectiveness of our DDM-DG method in solving such multi-scale problems.

Many multiscale methods are based on the idea of extracting macroscopic behavior of solutions by solving an array of microscale models over small domains. A key ingredient in such multiscale methods is the boundary condition and the size of the computational domain over which the microscale problems are solved. This problem is systematically investigated in the present paper in the context of modeling strongly heterogeneous media. Three different boundary conditions are considered: the periodic boundary condition, Dirichlet boundary condition, and the Neumann boundary condition. Each is applied to several benchmark problems: the random checker-board problem, periodic problem with isotropic macroscale behavior, periodic problem with anisotropic macroscale behavior and periodic laminated media. In each case, convergence studies are conducted as the domain size for the microscale problem is changed. Convergence rates as well as the size of fluctuations in the computed effective coefficients are compared for the different formulations. In addition, we will discuss a mixed Dirichlet-Neumann boundary condition that is often used in porous medium modeling. We explain why that leads to unsatisfactory results and how it can be corrected. Also discussed are the different averaging methods used in extracting the effective coefficients.

We present an efficient numerical algorithm for simulating chemical kinetic systems with multiple time scales. This algorithm is an improvement of the traditional stochastic simulation algorithm (SSA), also known as Gillespie's algorithm. It is in the form of a nested SSA and uses an outer SSA to simulate the slow reactions with rates computed from realizations of inner SSAs that simulate the fast reactions. The algorithm itself is quite general and seamless, and it amounts to a small modification of the original SSA. Our analysis of such multi-scale chemical kinetic systems allows us to identify the slow variables in the system, derive effective dynamics on the slow time scale, and provide error estimates for the nested SSA. Efficiency of the nested SSA is discussed using these error estimates, and illustrated through several numerical examples.

We discuss numerical methods for linear and nonlinear transport equations with multiscale velocity fields. These methods are themselves multiscaled in nature in the sense that they use macro and micro grids, multiscale test functions. We demonstrate the efficiency of these methods and apply them to two-phase flow in heterogeneous porous media.

The framework of the heterogeneous multiscale method (HMM) is used to develop numerical methods for the study of macroscale dynamics of fluids in situations, where either the constitutive relation or the boundary conditions are not explicitly available and have to be inferred from microscopic models such as molecular dynamics. Continuum hydrodynamics is used as the macroscopic model, while molecular dynamics serves as the microscopic model and is used to supply the necessary data, e.g., the stress or the boundary condition, for the macroscopic model. Scale separation is exploited so that the macroscopic variables can be evolved in macroscopic spatial/temporal scales using data that are estimated from molecular dynamics simulation on microscale spatial/temporal domains. This naturally decouples the micro and macrospatial and temporal scales whenever possible. Applications are presented for models of complex fluids, contact line dynamics, and a simple model of non-trivial fluid-solid interactions.

We present a hybrid method for simulating epitaxial growth that combines kinetic Monte-Carlo (KMC) simulations with the Burton-Cabrera-Frank model for crystal growth. This involves partitioning the computational domain into KMC regions and regions where we time-step a discretized diffusion equation. Computational speed and accuracy are discussed. We find that the method is significantly faster than KMC while accounting for stochastic fluctuations in a comparable way.

In current FFT-based algorithms for micromagnetics simulations, the boundary is typically replaced by a staircase approximation along the grid lines, either eliminating the incomplete cells or replacing them by complete cells. Sometimes the magnetizations at the boundary cells are weighted by the volume of the sample in the corresponding cell. We show that this leads to large errors in the computed exchange and stray fields. One consequence of this is that the predicted switching mechanism depends sensitively on the orientation of the numerical grid. We present a boundary-corrected algorithm to efficiently and accurately handle the incomplete cells at the boundary. We show that this boundary-corrected algorithm greatly improves the accuracy in micromagnetics simulations. We demonstrate by using A. Arrott's example of a hexagonal element that the switching mechanism is predicted independently of the grid orientation.

We present a coupled atomistic-continuum method for the modeling of defects and interface dynamics of crystalline materials. The method uses atomistic models such as molecular dynamics near defects and interfaces, and continuum models away from defects and interfaces. We propose a new class of matching conditions between the atomistic and continuum regions. These conditions ensure the accurate passage of large scale information between the atomistic and continuum regions and at the same time minimize the reflection of phonons at the atomistic-continuum interface. They can be made adaptive if we choose appropriate weight functions. We present applications to dislocation dynamics, friction between two-dimensional crystal surfaces and fracture dynamics. We compare results of the coupled method and the detailed atomistic model., Comment: 48 pages, 20 figures

One of the main difficulties in micromagnetics simulation is the severe time step constraint introduced by the exchange field. Using standard explicit integrators leads to a physical time step of sub-pico seconds, which is often two orders of magnitude smaller than the fastest physical time scales. Direct implicit integrators require solving complicated, coupled systems. In this paper, we introduce an implicit method whose complexity is comparable to solving the scalar heat equation implicitly. This method is based on a combination of a Gauss–Seidel implementation of a fractional step implicit solver for the gyromagnetic term, and the projection method for the heat flow of harmonic maps. This method allows us to carry out fully resolved calculations for the switching of the magnetization in micron-sized elements.

Simple, efficient, and accurate finite difference methods are introduced for 3D unsteady viscous incompressible flows in the vorticity?vector potential formulation on nonstaggered grids. Two different types of methods are discussed. They differ in the implementation of the normal component of the vorticity boundary condition and consequently the enforcement of the divergence free condition for vorticity. Both second-order and fourth-order accurate schemes are developed. A detailed accuracy test is performed, revealing the structure of the error and the effect of how the convective terms are discretized near the boundary. The influence of the divergence free condition for vorticity to the overall accuracy is studied. Results on the cubic driven cavity flow at Reynolds number 500 and 3200 are shown and compared with that of the MAC scheme.

We consider finite difference schemes based on the impulse density variable. We show that the original velocity?impulse density formulation of Oseledets is marginally ill-posed for the inviscid flow, and this has the consequence that some ordinarily stable numerical methods in other formulations become unstable in the velocity?impulse density formulation. We present numerical evidence of this instability. We then discuss the construction of stable finite difference schemes by requiring that at the numerical level the nonlinear terms be convertible to similar terms in the primitive variable formulation. Finally we give a simplified velocity?impulse density formulation which is free of these complications and yet retains the nice features of the original velocity?impulse density formulation with regard to the treatment of boundary. We present numerical results on this simplified formulation for the driven cavity flow on both the staggered and non-staggered grids.

A new fourth-order accurate finite difference scheme for the computation of unsteady viscous incompressible flows is introduced. The scheme is based on the vorticity-stream function formulation. It is essentially compact and has the nice features of a compact scheme with regard to the treatment of boundary conditions. It is also very efficient, at every time step or RungeÂ?Kutta stage, only two Poisson-like equations have to be solved. The Poisson-like equations are amenable to standard fast Poisson solvers usually designed for second order schemes. Detailed comparison with the second-order scheme shows the clear superiority of this new fourth-order scheme in resolving both the boundary layers and the gross features of the flow. This efficient fourth-order scheme also made it possible to compute the driven cavity flow at Reynolds number 106on a 10242grid at a reasonable cost. Fourth-order convergence is proved under mild regularity requirements. This is the first such result to our knowledge.

This paper discusses three basic issues related to the design of finite difference schemes forunsteadyviscous incompressible flows using vorticity formulations: the boundary condition for vorticity, an efficient time-stepping procedure, and the relation between these schemes and the ones based on velocity?pressure formulation. We show that many of the newly developed global vorticity boundary conditions can actually be written as some local formulas derived earlier. We also show that if we couple a standard centered difference scheme with third- or fourth-order explicit Runge?Kutta methods, the resulting schemes havenocell Reynolds number constraints. For high Reynolds number flows, these schemes are stable under the CFL condition given by the convective terms. Finally, we show that the classical MAC scheme is the same as Thom's formula coupled with second-order centered differences in the interior, in the sense that one can define discrete vorticity in a natural way for the MAC scheme and get the same values as the ones computed from Thom's formula. We use this to derive an efficient fourth-order Runge?Kutta time discretization for the MAC scheme from the one for Thom's formula. We present numerical results for driven cavity flow at high Reynolds number (105).