科研成果

2017
Liu W, Shao S, Li Z. Relativistic Explicit Correlation: Problems and Solutions. In: Handbook of Relativistic Quantum Chemistry. Berlin: Springer; 2017. pp. 531-545. 访问链接Abstract
The fundamental problems inherent in relativistic explicit correlation are highlighted, with practical suggestions for guiding future development of relativistic explicitly correlated wave function methods. 
Chang K-C, Shao S, Zhang D. Cheeger's cut, maxcut and the spectral theory of 1-Laplacian on graphs. SCIENCE CHINA Mathematics [Internet]. 2017;60(11):1963-1980. 访问链接Abstract
This is primarily an expository paper surveying up-to-date known results on the spectral theory of 1-Laplacian on graphs and its applications to the Cheeger cut, maxcut and multi-cut problems. The structure of eigenspace, nodal domains, multiplicities of eigenvalues, and algorithms for graph cuts are collected.
2016
Xiong Y, Chen Z, Shao S. An advective-spectral-mixed method for time-dependent many-body Wigner simulations. SIAM Journal on Scientific Computing [Internet]. 2016;38(4):B491–B520. 访问链接Abstract
As a phase space language for quantum mechanics, the Wigner function approach bears a close analogy to classical mechanics and has been drawing growing attention, especially in simulating quantum many-body systems. However, deterministic numerical solutions have been almost exclusively confined to one-dimensional one-body systems and few results are reported even for one-dimensional two-body problems. This paper serves as the first attempt to solve the time-dependent many-body Wigner equation through a grid-based advective-spectral-mixed method. The main feature of the method is to resolve the linear advection in $(\bm{x},t)$-space by an explicit three-step characteristic scheme coupled with the piecewise cubic spline interpolation, while the Chebyshev spectral element method in $\bm k$-space is adopted for accurate calculation of the nonlocal pseudo-differential term. Not only the time step of the resulting method is not restricted by the usual CFL condition and thus a large time step is allowed, but also the mass conservation can be maintained. In particular, for the system consisting of identical particles,  the advective-spectral-mixed method can also rigorously preserve physical symmetry relations. The performance is validated through several typical numerical experiments, like the Gaussian barrier scattering, electron-electron interaction and a Helium-like system, where the third-order accuracy against both grid spacing and time stepping is observed.  
Mertens FG, Cooper F, Quintero NR, Shao S, Khare A, Saxena A. Solitary waves in the nonlinear Dirac equation in the presence of external driving forces. Journal of Physics A: Mathematical and Theoretical [Internet]. 2016;49:065402. 访问链接Abstract
We consider the nonlinear Dirac (NLD) equation in (1 + 1) dimensions with scalar–scalar self interaction $\displaystyle \frac{{g}^{2}}{\kappa +1}{(\bar{{\rm{\Psi }}}{\rm{\Psi }})}^{\kappa +1}$ in the presence of external forces as well as damping of the form $f(x)-{\rm{i}}\mu {\gamma }^{0}{\rm{\Psi }}$, where both f and Ψ are two-component spinors. We develop an approximate variational approach using collective coordinates (CC) for studying the time dependent response of the solitary waves to these external forces. This approach predicts intrinsic oscillations of the solitary waves, i.e. the amplitude, width and phase all oscillate with the same frequency. The translational motion is also affected, because the soliton position oscillates around a mean trajectory. For $\kappa =1$ we solve explicitly the CC equations of the variational approximation for slow moving solitary waves in a constant external force without damping and find reasonable agreement with solving numerically the CC equations. We then compare the results of the variational approximation with no damping with numerical simulations of the NLD equation for $\kappa =1$, when the components of the external force are of the form ${f}_{j}={r}_{j}\mathrm{exp}(-{\rm{i}}{Kx})$ and again find agreement if we take into account a certain linear excitation with specific wavenumber that is excited together with the intrinsic oscillations such that the momentum in a transformed NLD equation is conserved.
Wang D, Zhou Y, Shao S. Efficient implementation of smoothed particle hydrodynamics (SPH) with plane sweep algorithm. Communications in Computational Physics [Internet]. 2016;19(3):770-800. 访问链接Abstract
Neighbour search (NS) is the core of any implementations of smoothed particle hydrodynamics (SPH). In this paper, we present an efficient $\mathcal{O}(N\log N)$ neighbour search method based on the plane sweep (PW) algorithm with $N$ being the number of SPH particles. The resulting method, dubbed the PWNS method, is totally independent of grids (i.e., purely meshfree) and capable of treating variable smoothing length, arbitrary particle distribution and heterogenous kernels. Several state-of-the-art data structures and algorithms, e.g., the segment tree and the Morton code, are optimized and implemented. By simply allowing multiple lines to sweep the SPH particles simultaneously from different initial positions, a parallelization of the PWNS method with satisfactory speedup and load-balancing can be easily achieved. That is, the PWNS SPH solver has a great potential for large scale fluid dynamics simulations.
2015
Shao S, Sellier JM. Comparison of deterministic and stochastic methods for time-dependent Wigner simulations. Journal of Computational Physics [Internet]. 2015;300:167-185. 访问链接Abstract
Recently a Monte Carlo method based on signed particles for time-dependent simulations of the Wigner equation has been proposed. While it has been thoroughly validated against physical benchmarks, no technical study about its numerical accuracy has been performed. To this end, this paper presents the first step towards the construction of firm mathematical foundations for the signed particle Wigner Monte Carlo method. An initial investigation is performed by means of comparisons with a cell average spectral element method, which is a highly accurate deterministic method and utilized to provide reference solutions. Several different numerical tests involving the time-dependent evolution of a quantum wave-packet are performed and discussed in deep details. In particular, this allows us to depict a set of crucial criteria for the signed particle Wigner Monte Carlo method to achieve a satisfactory accuracy.  
Xu J, Shao S, Tang H, Wei D. Multi-hump solitary waves of a nonlinear Dirac equation. Communications in Mathematical Sciences [Internet]. 2015;13(5):1219-1242. 访问链接Abstract
This paper concentrates on a (1+1)-dimensional nonlinear Dirac (NLD) equation with a general self-interaction, being a linear combination of the scalar, pseudoscalar, vector and axial vector self-interactions to the power of the integer $k+1$. The solitary wave solutions to the NLD equation are analytically derived, and the upper bounds of the hump number in the charge, energy and momentum densities for the solitary waves are proved analytically in theory. The results show that: (1) for a given integer $k$, the hump number in the charge density is not bigger than $4$, while that in the energy density is not bigger than $3$; (2) those upper bounds can only be achieved in the situation of higher nonlinearity, namely, $k\in\{5,6,7,\cdots \}$ for the charge density and $k\in\{3,5,7,\cdots\}$ for the energy density; (3) the momentum density has the same multi-hump structure as the energy density; (4) more than two humps (resp. one hump) in the charge (resp. energy) density can only happen under the linear combination of the pseudoscalar self-interaction and at least one of the scalar and vector (or axial vector) self-interactions. Our results on the multi-hump structure will be interesting in the interaction dynamics for the NLD solitary waves.
Chang K-C, Shao S, Zhang D. The 1-Laplacian Cheeger cut: Theory and algorithms. Journal of Computational Mathematics [Internet]. 2015;33(5):443-467. 访问链接Abstract
This paper presents a detailed review of both theory and algorithms for the Cheeger cut based on the graph $1$-Laplacian. In virtue of the cell structure of the feasible set, we propose a cell descend (CD) framework for achieving the Cheeger cut. While plugging the relaxation to guarantee the decrease of the objective value in the feasible set, from which both the inverse power (IP) method and the steepest descent (SD) method can also be recovered, we are able to get two specified CD methods. Comparisons of all these methods are conducted on several typical graphs.
2014
Shao S, Quintero NR, Mertens FG, Cooper F, Khare A, Saxena A. Stability of solitary waves in the nonlinear Dirac equation with arbitrary nonlinearity. Physical Review E [Internet]. 2014;90:032915. 访问链接Abstract
We consider the nonlinear Dirac equation in 1+1 dimension with scalar-scalar  self interaction $ \frac{ g^2}{ \kappa+1} ( {\bar \Psi} \Psi)^{ \kappa+1}$ and with mass $m$.    Using  the exact analytic form for  rest frame solitary waves of the form $\Psi(x,t) = \psi(x) e^{-i \omega t}$  for arbitrary $ \kappa$, we  discuss the validity  of various approaches to understanding stability that were successful for the nonlinear Schr\"odinger equation. In particular we study the validity of a version of  Derrick's theorem, the criterion of Bogolubsky  as well as  the Vakhitov-Kolokolov criterion, and find that these criteria yield inconsistent results. Therefore, we study the stability by numerical  simulations using a recently developed 4th-order operator splitting integration method. For different ranges of $\kappa$ we map out the stability regimes in $\omega$. We find that all stable nonlinear Dirac solitary waves have a one-hump profile, but not all one-hump waves are stable, while all waves with two humps are unstable. We also find that the time $t_c$, it takes for the instability to set in, is an exponentially increasing function of $\omega$ and $t_c$ decreases monotonically with increasing $\kappa$.
Lin L, Lu J, Shao S. Analysis of the time reversible Born-Oppenheimer molecular dynamics. Entropy [Internet]. 2014;16:110-137. 访问链接Abstract
We analyze the time reversible Born-Oppenheimer molecular dynamics (TRBOMD) scheme, which preserves the time reversibility of the Born-Oppenheimer molecular dynamics even with non-convergent self-consistent field iteration. In the linear response regime, we derive the stability condition, as well as the accuracy of TRBOMD for computing physical properties, such as the phonon frequency obtained from the molecular dynamics simulation. We connect and compare TRBOMD with Car-Parrinello molecular dynamics in terms of accuracy and stability. We further discuss the accuracy of TRBOMD beyond the linear response regime for non-equilibrium dynamics of nuclei. Our results are demonstrated through numerical experiments using a simplified one-dimensional model for Kohn-Sham density functional theory. 
Wang D, Shao S, Yan C, Cai W, Zeng X. Feature-scale simulations of particulate slurry flows in chemical mechanical polishing by smoothed particle hydrodynamics. Communications in Computational Physics [Internet]. 2014;16(5):1389-1418. 访问链接Abstract
In this paper, the mechanisms of material removal in chemical mechanical polishing (CMP) processes are investigated in detail by the smoothed particle hydrodynamics (SPH) method. The feature-scale behaviours of slurry flow, rough pad, wafer defects, moving solid boundaries, slurry-abrasive interactions, and abrasive collisions are modelled and simulated. Compared with previous work on CMP simulations, our simulations incorporate more realistic physical aspects of the CMP process, especially the effect of abrasive concentration in the slurry flows. The preliminary results on slurry flow in CMP provide microscopic insights on the experimental data of the relation between the removal rate and abrasive concentration and demonstrate that SPH is a suitable method for the research of CMP processes.
2013
Xu J, Shao S, Tang H. Numerical methods for nonlinear Dirac equation. Journal of Computational Physics [Internet]. 2013;245:131-149. 访问链接Abstract
This paper presents a review of the current state-of-the-art of numerical methods for nonlinear Dirac (NLD) equation. Several methods are extendedly proposed for the (1+1)-dimensional NLD equation with the scalar and vector self-interaction and analyzed in the way of the accuracy and the time reversibility as well as the conservation of the discrete charge, energy and linear momentum. Those methods are the Crank-Nicolson (CN) schemes, the linearized CN schemes, the odd-even hopscotch scheme, the leapfrog scheme, a semi-implicit finite difference scheme, and the exponential operator splitting (OS) schemes. The nonlinear subproblems resulted from the OS schemes are analytically solved by fully exploiting the local conservation laws of the NLD equation. The effectiveness of the various numerical methods, with special focus on the error growth and the computational cost, is illustrated on two numerical experiments, compared to two high-order accurate Runge-Kutta discontinuous Galerkin methods. Theoretical and numerical comparisons show that the high-order accurate OS schemes may compete well with other numerical schemes discussed here in terms of the accuracy and the efficiency. A fourth-order accurate OS scheme is further applied to investigating the interaction dynamics of the NLD solitary waves under the scalar and vector self-interaction. The results show that the interaction dynamics of two NLD solitary waves   depend on the exponent power of the self-interaction in the NLD equation; collapse happens after collision of two equal one-humped NLD solitary waves under the cubic vector self-interaction in contrast to no collapse scattering for corresponding quadric case. 
Lin L, Shao S, E W. Efficient iterative method for solving the Dirac-Kohn-Sham density functional theory. Journal of Computational Physics [Internet]. 2013;245:205-217. 访问链接Abstract
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. 
2012
Li Z, Shao S, Liu W. Relativistic explicit correlation: Coalescence conditions and practical suggestions. Journal of Chemical Physics [Internet]. 2012;136:144117. 访问链接Abstract
To set up the general framework for relativistic explicitly correlated wave function methods, the electron-electron coalescence conditions are derived for the wave functions of the Dirac-Coulomb (DC), Dirac-Coulomb-Gaunt (DCG), Dirac-Coulomb-Breit (DCB), modified Dirac-Coulomb (MDC), and zeroth-order regularly approximated (ZORA) Hamiltonians. The manipulations make full use of the internal symmetries of the reduced two-electron Hamiltonians such that the asymptotic behaviors of the wave functions emerge naturally. The results show that, at the coalescence point of two electrons, the wave functions of the DCG Hamiltonian are regular, while those of the DC and DCB Hamiltonians have weak singularities of the type $r_{12}^{\nu}$ with $\nu$ being negative and of $\mathcal{O}(\alpha^2)$. The behaviors of the MDC wave functions are related to the original ones in a simple manner, while the spin-free counterparts are somewhat different due to the complicated electron-electron interaction. The behaviors of the ZORA wave functions depend on the chosen potential in the kinetic energy operator. In the case of the nuclear attraction, the behaviors of the ZORA wave functions are very similar to those of the nonrelativistic ones, just with an additional correction of $\mathcal{O}(\alpha^2)$ to the nonrelativistic cusp condition. However, if the Coulomb interaction is also included, the ZORA wave functions become close to the large-large components of the DC wave functions. Note that such asymptotic expansions of the relativistic wave functions are only valid within an extremely small convergence radius $R_c$ of $\mathcal{O}(\alpha^2)$. Beyond this radius, the behaviors of the relativistic wave functions are still dominated by the nonrelativistic limit, as can be seen in terms of direct perturbation theory (DPT) of relativity. However, as the two limits $\alpha\rightarrow0$ and $r_{12}\rightarrow0$ do not commute, DPT is doomed to fail due to incorrect descriptions of the small-small component $\Psi^{SS}$ of the DC wave function for $r_{12}<R_c$. Another deduction from the possible divergence of  $\Psi^{SS}$ at $r_{12}=R_c$ is that the DC Hamiltonian has no bound electronic states, although the last word cannot be said. These findings enrich our understandings of relativistic wave functions. On the practical side, it is shown that, under the no-pair approximation, relativistic explicitly correlated wave function methods can be made completely parallel to the nonrelativistic counterparts, as demonstrated explicitly for MP2-F12. Yet, this can only be achieved by using an extended no-pair projector.
Shao S, Qian T. A variational model for two-phase immiscible electro-osmotic flows. Communications in Computational Physics [Internet]. 2012;11(3):831-862. 访问链接Abstract
We develop a continuum hydrodynamic model for two-phase immiscible flows that involve electroosmotic effect in an electrolyte and moving contact line at solid surfaces. The model is derived through a variational approach based on the On- sager principle of minimum energy dissipation. This approach was first presented in the derivation of a continuum hydrodynamic model for moving contact line in neu- tral two-phase immiscible flows (Qian, Wang, and Sheng, J. Fluid Mech. 564, 333–360 (2006)). Physically, the electroosmotic effect can be formulated by the Onsager prin- ciple as well in the linear response regime. Therefore, the same variational approach is applied here to the derivation of the continuum hydrodynamic model for charged two-phase immiscible flows where one fluid component is an electrolyte exhibiting electroosmotic effect on a charged surface. A phase field is employed to model the diffuse interface between two immiscible fluid components, one being the electrolyte and the other a nonconductive fluid, both allowed to slip at solid surfaces. Our model consists of the incompressible Navier-Stokes equation for momentum transport, the Nernst-Planck equation for ion transport, the Cahn-Hilliard phase-field equation for interface motion, and the Poisson equation for electric potential, along with all the necessary boundary conditions. In particular, all the dynamic boundary conditions at solid surfaces, including the generalized Navier boundary condition for slip, are de- rived together with the equations of motion in the bulk region. Numerical examples in two-dimensional space, which involve overlapped electric double layer fields, have been presented to demonstrate the validity and applicability of the model, and a few salient features of the two-phase immiscible electroosmotic flows at solid surface. The wall slip in the vicinity of moving contact line and the Smoluchowski slip in the electric double layer are both investigated.
2011
Shao S, Lu T, Cai W. Adaptive conservative cell average spectral element methods for transient Wigner equation in quantum transport. Communications in Computational Physics [Internet]. 2011;9(3):711-739. 访问链接Abstract
A new adaptive cell average spectral element method (SEM) is proposed to solve the time-dependent Wigner equation for transport in quantum devices. The proposed cell average SEM allows adaptive non-uniform meshes in phase spaces to reduce the high-dimensional computational cost of Wigner functions while preserving exactly the mass conservation for the numerical solutions. The key feature of the pro- posed method is an analytical relation between the cell averages of the Wigner function in the k-space (local electron density for finite range velocity) and the point values of the distribution, resulting in fast transforms between the local electron density and local fluxes of the discretized Wigner equation via the fast sine and cosine transforms. Numerical results with the proposed method are provided to demonstrate its high accuracy, conservation, convergence and a reduction of the cost using adaptive meshes.
2008
Jiang H, Shao S, Cai W, Zhang P. Boundary treatments in non-equilibrium Green's function (NEGF) methods for quantum transport in nano-MOSFETs. Journal of Computational Physics [Internet]. 2008;227:6553–6573. 访问链接Abstract
Non-equilibrium Green’s function (NEGF) is a general method for modeling non-equilibrium quantum transport in open mesoscopic systems with many body scattering effects. In this paper, we present a unified treatment of quantum device boundaries in the framework of NEGF with both finite difference and finite element discretizations. Boundary treat- ments for both types of numerical methods, and the resulting self-energy R for the NEGF formulism, representing the dis- sipative effects of device contacts on the transport, are derived using auxiliary Green’s functions for the exterior of the quantum devices. Numerical results with both discretization schemes for an one-dimensional nano-device and a 29 nm double gated MOSFET are provided to demonstrate the accuracy and flexibility of the proposed boundary treatments.
Cai W, Ji X, Sun J, Shao S. A Schwarz generalized eigen-oscillation spectral element method (GeSEM) for 2-D high frequency electromagnetic scattering in dispersive inhomogeneous media. Journal of Computational Physics [Internet]. 2008;227(23):9933-9954. 访问链接Abstract
In this paper, we propose a parallel Schwarz generalized eigen-oscillation spectral element method (GeSEM) for 2-D complex Helmholtz equations in high frequency wave scattering in dispersive inhomogeneous media. This method is based on the spectral expansion of complex generalized eigen-oscillations for the electromagnetic fields and the Schwarz non-overlapping domain decomposition iteration method. The GeSEM takes advantages of a special real orthogonality property of the complex eigen-oscillations and a new radiation interface condition for the system of equations for the spectral expansion coefficients. Numerical results validate the high resolution and the flexibility of the method for various materials.
Shao S, Tang H. Interaction of solitary waves with a phase shift in a nonlinear Dirac model. Communications in Computational Physics [Internet]. 2008;3(4):950-967. 访问链接Abstract
This paper presents a further numerical study of the interaction dynamics for solitary waves in a nonlinear Dirac model with scalar self-interaction, the Soler model, by using a fourth order accurate Runge-Kutta discontinuous Galerkin method. The phase plane method is employed for the first time to analyze the interaction of Dirac solitary waves and reveals that the relative phase of those waves may vary with the interaction. In general, the interaction of Dirac solitary waves depends on the initial phase shift. If two equal solitary waves are in-phase or out-of-phase initially, so are they during the interaction; if the initial phase shift is far away from 0 and π, the relative phase begins to periodically evolve after a finite time. In the interaction of out-of-phase Dirac solitary waves, we can observe: (a) full repulsion in binary and ternary collisions, depending on the distance between initial waves; (b) repulsing first, attracting afterwards, and then collapse in binary and ternary collisions of initially resting two-humped waves; (c) one-overlap interaction and two-overlap interaction in ternary collisions of initially resting waves.
2006
Shao S, Cai W, Tang H. Accurate calculation of Green's function of the Schrödinger equation in a block layered potential. Journal of Computational Physics [Internet]. 2006;219(2):733-748. 访问链接Abstract
In this paper a new algorithm is presented for calculating the Green’s function of the Schrödinger equation in the presence of block layered potentials. Such Green’s functions have various and practical applications in quantum modelling of electron transport within nano-MOSFET transistors. The proposed method is based on expansions of the eigenfunctions of the subordinate Sturm–Liouville problems and a collocation matching procedure along possibly curved interfaces of the potential blocks. Accurate numerical results are provided to validate the proposed algorithm.

Pages