\+<EMAIL_ADDRESS>
# Implementation of LDA+DMFT with pseudo-potential-plane-wave method

Jian-Zhou Zhao1,2, Jia-Ning Zhuang1, Xiao-Yu Deng1, Yan Bi2, Ling-Cang Cai2,
Zhong Fang1, and Xi Dai1 1Institute of Physics, Chinese Academy of Sciences,
Beijing 100190, People’s Republic of China  2National key Laboratory of Shock
Wave and Detonation Physics, Institute of Fluid Physics, China Academy of
Engineering Physics, Mianyang, 621900, China

###### Abstract

In this paper, we propose an efficient implementation of combining Dynamical
Mean field theory (DMFT) with electronic structure calculation based on the
local density approximation (LDA). The pseudo-potential-plane-wave method is
used in the LDA part, which makes it possible to be applied to large systems.
The full loop self consistency of the charge density has been reached in our
implementation which allows us to compute the total energy related properties.
The procedure of LDA+DMFT is introduced in detail with a complete flow chart.
We have also applied our code to study the electronic structure of several
typical strong correlated materials, including Cerium, Americium and NiO. Our
results fit quite well with both the experimental data and previous studies.

###### pacs:

71.27.+a, 71.15.Mb, 71.15.Nc, 71.20.-b

##  1 Introduction

The first principle calculation based on the density functional theory (DFT)
with the local density approximation (LDA) and its generalization generalized
gradient approximation (GGA) is very successful in predicting ground-state
properties and band structures of a wide range of real materials. However it
is known for a long time that it is not sufficient to calculate the electronic
structure of those strongly correlated materials (i.e. transition metal
oxides, actinide and lanthanide-based materials, high $T_{c}$ superconductors)
by these methods alone even on the qualitative level. In order to overcome
these shortcomings of the traditional DFT-LDA scheme, remedies such as LDA+$U$
have been proposed, which can describe some of the strong correlated materials
with long range order in their ground states. While even LDA+$U$ can not be
applied to many of the strongly correlated materials, for example, the
paramagnetic Mott insulator phase and the materials containing unfilled
f-shell with strong multiplet effects.

Therefore it is important to have a first principle method which can be
applied to the strongly correlated materials with the ability to capture the
dynamical properties, finite temperature properties and multiplet effects. In
the last two decades, the dynamical mean field theory (DMFT) has been
developed very fast to be the standard tool to deal with the on-site
correlation effect in the limit of large dimension [1]. After being
successfully applied to many model systems, DMFT has been considered as a
powerful tool to capture the on-site correlation effects based on the Hubbard
like Hamiltonians containing both the local interaction terms and the single
particle Hamiltonians extracted from LDA. Therefore the LDA+DMFT method which
combines the DFT-based band structure techniques with DMFT has been proposed
and developed quickly in the last decade. By DMFT the local correlation effect
can be well described by the self energy, which has frequency dependence and
in general takes the matrix form within the subspace spanned by the
correlation orbitals. Using the Green’s function containing self energy, many
of the physical properties of strongly correlated materials can be calculated,
i.e. the electronic spectral function, the total energy, the optical
conductivity and the local spin susceptibility.

Most of the LDA+DMFT calculation till now have been performed by the partially
self-consistent scheme, where the local self-energy is obtained by the DMFT
calculation with a fixed LDA Hamiltonian generated from the fixed LDA charge
density. Therefore, in this simplified scheme one neglects the inference of
the strong on-site Coulomb interaction on the charge density. The above
mentioned “one-shot” DMFT calculation works quite well for the electronic
structure. While in order to obtain reliable total energy related properties,
i.e. the equilibrium volume, the elastic constants and the phonon frequencies,
the LDA+DMFT scheme with full charge density self consistency is needed. Up to
date, there have been several fully self-consistent LDA+DMFT schemes [2, 3] as
well as actual implementations [4, 5, 6, 7]. As mentioned in [7], there are
three major issues that have to be addressed in DFT+DMFT implementations: i)
quality of the basis set, ii) quality of the impurity solvers, iii) choice of
correlated orbitals onto which the full Green’s function is projected. In the
implementations of LDA+DMFT mentioned above with full charge self-consistency
, the basis set of the linear muffin-tin orbitals (LMTO) are used, and the
local correlated orbitals are naturally chosen based on the muffin-tin
orbitals . While on the other hand, there are very few implementations of
LDA+DMFT with full charge self-consistency are based on the plane wave
methods.

In the present paper, we implement a full self-consistent LDA+DMFT scheme
using pseudo-potential-plane-wave as the basis set, and we use either atomic
orbitals or Wannier orbitals depending on different physics systems. The
Hubbard-I approximation is used as the impurity solver in the present paper,
but it can be replace by more precise solver like continuous time quantum
Monte Carlo (CTQMC) for metallic systems. The present paper is organized as
follows. In section 2, we describe the way to choose and construct the
correlated orbitals onto which the local interactions exert. In section 3, we
show the full-loop flow and LDA+DMFT(Hubbard-I) formalism in detail. We apply
this LDA+DMFT approach to several correlated materials such as Ce, Am, and NiO
in section 4. Finally we conclude our work in section 5.

##  2 Projection onto localized orbitals

In the LDA+DMFT calculation, all the energy bands are divided into two groups.
And only the on-site interactions among those local orbitals are treated more
precisely by DMFT. Therefore the choice of localized orbitals does affect the
results obtained by LDA+DMFT method and becomes one of the important issues in
LDA+DMFT. A natural choice for the local orbitals is a set of atomic like wave
functions with corresponding $d$ or $f$ characters. Typical atomic like local
orbitals is the linear muffin-tin orbitals (LMTOs)[8] adopted in early
LDA+DMFT implementations[9, 10, 11]. However a more physical choice should be
wannier functions, since the shape of the local orbitals will be altered in
crystals especially when there is strong hybridization between localized
orbital and delocalized $p$-type or $s$-type orbitals. The wannier functions
are not uniquely determined by the Bloch wave functions, so different choice
can be made, for example, the Nth-order muffin-tin orbitals[12], the projected
wannier functions[13] and the Maximally localized Wannier functions
(MLWFs)[14, 15]. Comparisons of these choices have been made in previous
literature. In this paper, two kinds of local orbitals, atomic like functions
and the projected wannier functions are adopted according to the different
systems.

This implementation of LDA+DMFT method reported in this paper is developed on
an existing DFT package BSTATE (Beijing Simulation Tool of Atomic TEchnique),
which is based on the pseudo-potentials and plane waves. Unlike the LMTO
methods, the local orbitals do not enter the basis set of plane waves and
should be constructed in a suitable way. Projected wannier functions or MLWFs
as local orbitals have been used by previous reported LDA+DMFT implementations
based on pseudo-potentials (or projector augmented wave (PAW)) method) and
plane wave scheme[16, 17, 18]. In this report, two kinds of orbitals are used
according to the need, one is directly derived from the atomic wave function
of an isolated atom and the other is projected wannier functions constructed
from these atomic wave functions. In order to be self-contained the
construction procedure is presented below in detail. Since all our
calculations are based on plane wave set, the local orbitals will be projected
onto these plane waves.

First we consider the atomic like local orbitals. The localized nature of the
correlated bands with $d$ or $f$ characters in crystals ensures that these
local orbitals are very similar to the corresponding atomic wave functions of
an isolated atom. The atomic wave functions can be picked as the local
orbitals directly if the hybridization in crystal could be neglected. All
atomic wave functions for the isolated atom could be obtained by the solving
the all-electron radical Schrödinger equation

$(T+V(r))\psi_{nl}^{ALL}(r)=E_{nl}\psi_{nl}^{ALL}(r)$ (1)

In this manner, all orbitals are well defined by the primer quantum number $n$
and angular quantum number $l$. In realistic systems, the local orbitals
usually come from the partially filled $d$ or $f$ shell and could be picked
according to its character. Often there is only one type of local orbitals, so
the local subspace can be labeled by $l$. Of course, considering the spherical
symmetry of the isolated atom, the spherical harmonics should be multiplied
the radical wave functions to form a complete set. The local orbitals are

$|\phi_{lm,I}\rangle=|\psi_{l,I}^{ALL}Y_{lm}\rangle,$ (2)

and in which $m=1\dots 2l+1$ denote different angular components.

In methods based on plane waves, it is not desirable to use the all-electron
wave function directly since it requires a large amount of plane waves to
expand in momentum space. To avoid this problem, in PAW method all-electron
atomic partial waves can be used, while in pseudo-potentials method, pseudo
atomic wave functions can be chosen. The latter is used in this paper. During
the generation of pseudo-potentials, the Schrödinger equation obeyed by the
pseudo wave functions is as below

$(T+V^{PS}(r))\psi_{nl}^{PS}(r)=E_{nl}^{PS}\psi_{nl}^{PS}(r)$ (3)

The details on the generation of pseudo-potentials can be referred to previous
literature[19, 20]. The spirit of pseudo-potentials is quite straight forward.
Beyond a core radius $r_{c}$, the pseudo-potential $V^{PS}$ play the role of a
scattering center just as a real atomic potential $V^{ALL}$ do. Thus, the
pseudo eigen energy $E_{nl}^{PS}$ should be the same as the realistic one
$E_{nl}^{ALL}$, and both the pseudo wave functions $|\psi_{n,I}^{PS}\rangle$
and the pseudo-potential $V^{PS}$ coincide with the exact wave functions
$|\psi_{n,I}^{ALL}\rangle$ and the realistic potential $V^{ALL}$ of an
isolated atoms beyond this core radius $r_{c}$, respectively.

$\eqalign{E_{nl}^{PS}=E_{nl};\cr
V^{PS}(r)=V(r),\qquad\psi^{PS}(r)=\psi^{ALL}(r)\qquad r>r_{c}}$ (4)

The quality of pseudo wave functions is the same as the quality of the pseudo-
potentials, which could be justified by comparing simple DFT calculations with
accepted results. The pseudo wave functions bear the same atomic features as
the exact atomic wave functions, and can be picked as the local orbitals as
above (angular quantum number $l$ ignored since usually only one kind of local
orbital are considered).

$|\phi_{m,I}\rangle=|\psi_{l,I}^{PS}Y_{lm}\rangle,$ (5)

It is convenient to transform the local orbitals from real space to momentum
space in crystal calculations

$|\phi_{m,{\mathbf{k}}}\rangle=\sum_{I}e^{i{\mathbf{k}}\cdot{\mathbf{R}}_{I}}|\phi_{m,I}\rangle.$
(6)

Since the atomic wave functions at different site $I$ are not orthogonal, an
orthogonalization procedure is needed for the above local orbitals.

The DFT calculations give a set of Bloch waves $|\Psi_{n\mathbf{k}}\rangle$
spanning the total Hilbert space in which lies the local subspace spanning by
the local orbitals. Thus, the physical choice of local orbitals is Wannier
functions derived from the Bloch waves. The Wannier functions are localized in
real space and could be constructed from the Bloch waves via a unitary
transformation. However, this unitary transformation is not unique since the
phase factors of Bloch waves are uncertain, which results in that there are
many kinds of Wannier functions in use, e.g., the projected Wannier functions,
NMTOs and MLWFs. Also it is not necessary to include all the Bloch waves but
the relevant bands lying in a certain energy range to construct the Wannier
functions. The local orbitals could be picked as a subset of these Wannier
functions with specified $d$ or $f$ localized character.

The projecting Wannier functions method is a simple way to construct wannier
functions, in which trial wave functions are projecting onto physical relevant
Bloch bands. The atomic like local orbitals above are used as trial functions
here. The Bloch bands can be selected by band indices ($N_{i},...N_{j}$) or by
an energy interval ($E_{i},E_{j}$) enclosing them.

$\displaystyle|W_{m,{\mathbf{k}}}\rangle$ $\displaystyle=$
$\displaystyle\sum_{n=N_{i}}^{N_{j}}|\psi_{n{\mathbf{k}}}\rangle\langle\psi_{n{\mathbf{k}}}|\phi_{lm,{\mathbf{k}}}\rangle$
(7) $\displaystyle=$
$\displaystyle\sum_{n(E_{i}<\epsilon_{n{\mathbf{k}}}<E_{j})}|\psi_{n{\mathbf{k}}}\rangle\langle\psi_{n{\mathbf{k}}}|\phi_{m,{\mathbf{k}}}\rangle.$

These orbitals need normalization and orthogonalization to be true Wannier
functions.

$|\tilde{W}_{m,{\mathbf{k}}}\rangle=\sum_{m^{\prime}}O_{m,m^{\prime}}({\mathbf{k}})^{-1/2}|{W}_{m,{\mathbf{k}}}\rangle$
(8)

in which $O$ is the overlap matrix between different
$|W_{m,{\mathbf{k}}}\rangle$s

$O_{m,m^{\prime}}({\mathbf{k}})=\langle
W_{m,{\mathbf{k}}}|W_{m^{\prime},{\mathbf{k}}}\rangle$ (9)

Although the construction method can not give the most localized Wannier
orbitals, the basic features of the localized bands are captured to a very
good extend as indicated by a few reports[13, 16, 17, 18], and also proved by
our calculations in this paper.

##  3 LDA+DMFT(Hubbard-I) Formalism

We will introduce the detail procedge of LDA+DMFT with full loop charge
density self consistency in this section. First we plot the general flow chart
in figure 1, which contains a inner DMFT loop and a outer charge density loop.
However, if we use Hubbard-I approximation as the impurity solver, the inner
DMFT loop is thus neglected.

Figure 1: The most general flow chart for LDA+DMFT scheme. The first LDA
calculation gives out the band structure on LDA level, and constructs local
orbitals, then the full Hamiltonian is established and solved by a DMFT loop,
after which the charge density can be recalculated by Fourier transforming the
k-space density matrix to real space. The “non-interacting” Hamiltonian is
thus regenerated based on the new charge density profile and the calculation
will be completed when the self consistency has been reached for the full
loop.

In the following subsections, we are going to describe the whole process in
detail.

###  3.1 First LDA Calculation

The first step of the full loop LDA+DMFT is an LDA self-consistent
calculation, whose main purpose is to generate an effective single particle
Hamiltonian $\hat{H}_{LDA}$ and construct the correlated orbitals.

Generally the effective LDA Hamiltonian can be expressed as the following:

$\displaystyle\hat{H}_{LDA}=\sum_{n\mathbf{k}}E_{n\mathbf{k}}\hat{c}^{\dagger}_{n\mathbf{k}}\hat{c}_{n\mathbf{k}}$
(10)

In the above equation, $n=1\sim N_{band}$ are the joint indices of band and
spin; $E_{n\mathbf{k}}$ is the eigen energy of the Kohn-Sham equation
determined by
$\hat{H}_{LDA}|\varphi_{n\mathbf{k}}\rangle=E_{n\mathbf{k}}|\varphi_{n\mathbf{k}}\rangle$,
where $|\varphi_{n\mathbf{k}}\rangle$ represents a set of orthonormal Bloch
functions.

The LDA calculation also complement the construction of local correlated
orbitals. In the present implementation the atomic orbitals and the Wannier
functions are two types of commonly used local basis. In LDA+DMFT, the local
interactions have been considered on the DMFT level only within the local
orbitals and in order to set up the DMFT self consistent equation we need to
obtain the overlap matrix between local orbitals and Bloch wave functions,
which takes the form of
$S^{\mathbf{k}}_{\alpha,n}=\langle\alpha_{\mathbf{k}}|\varphi_{n\mathbf{k}}\rangle$,
with Greek letter $\alpha$ denote the index of local orbitals. The
completeness of the Bloch basis set gives

$\displaystyle\sum_{n}|S^{\mathbf{k}}_{\alpha,n}|^{2}=1$ (11)

for any local orbital $\alpha$ and any k-point.

###  3.2 DMFT Loop

The purpose of the DMFT loop is to calculate the self-energy caused by the
local interactions through the DMFT self consistent loop. As we have mentioned
before, the DMFT loop is not necessity if we use Hubbard-I approximation as
the impurity solver. However, here we first introduce the algorithm with more
general impurity solvers.

By adding the local interaction terms to LDA, we get the total Hamiltonian of
the system which is to be solved by DMFT,

$\displaystyle\hat{H}_{LDA+DMFT}=\hat{H}_{LDA}+\sum_{i}(\hat{H}^{i}_{U}-\hat{V}^{i}_{DC})-\mu\hat{N}$
(12)

where $i$ is the site index. For a specified site, we remove the $i$ index and

$\displaystyle\hat{H}_{U}=\frac{1}{2}\sum_{\alpha\beta\gamma\delta}U_{\alpha\beta\gamma\delta}\hat{f}^{\dagger}_{\alpha}\hat{f}^{\dagger}_{\beta}\hat{f}_{\gamma}\hat{f}_{\delta}$
(13)

is written as a general form of two-body local interactions, in which the
creation and annihilation operators of correlated orbitals
$\hat{f}^{\dagger}_{\alpha}$ and $\hat{f}_{\alpha}$ are associated with Bloch
operators in term of

$\displaystyle\hat{f}^{\dagger}_{\alpha}=\sum_{n\mathbf{k}}\hat{c}^{\dagger}_{n\mathbf{k}}\langle\varphi_{n\mathbf{k}}|\alpha_{\mathbf{k}}\rangle$
(14)
$\displaystyle\hat{f}_{\alpha}=\sum_{n\mathbf{k}}\hat{c}_{n\mathbf{k}}\langle\alpha_{\mathbf{k}}|\varphi_{n\mathbf{k}}\rangle$
(15)

The third term of (12) is the double-counting term correspond to the
correlated energy that has already been considered in LDA calculation at a
Hartree-Fock mean field level,

$\hat{V}_{DC}=\sum_{\alpha}V_{DC}^{\alpha}\hat{f}^{\dagger}_{\alpha}\hat{f}_{\alpha}$
(16)

and this term will be discussed in section 3.2.3. The last term of (12)

$\hat{N}=\sum_{n\mathbf{k}}\hat{c}^{\dagger}_{n\mathbf{k}}\hat{c}_{n\mathbf{k}}$
(17)

is the total number of particle operator, while the chemical potential $\mu$
controls the occupation number of the unit cell.

####  3.2.1 Quantum Impurity Hamiltonian

In DMFT, the correlation problem on the lattice can be mapped to a quantum
impurity model, which contains the same on-site interaction and reads,

$\displaystyle\hat{H}_{imp}=\sum_{q\alpha}\epsilon_{q\alpha}\hat{c}^{\dagger}_{q\alpha}\hat{c}_{q\alpha}+\sum_{q\alpha}V_{q\alpha}(\hat{c}^{\dagger}_{q\alpha}\hat{f}_{\alpha}+h.c.)+\sum_{\alpha}E^{imp}_{\alpha}\hat{f}^{\dagger}_{\alpha}\hat{f}_{\alpha}+\hat{H}_{U}$
(18)

The inference from the rest of the lattice site besides the one considered in
the impurity model is simulated by a non-interacting “heat bath”, which is
described by the first term in (18). The second term describes the coupling
between the impurity site and the heat bath and the rest two terms describe
the local interactions.

####  3.2.2 Hybridization function and Weiss field

The above quantum impurity model can be solved by the impurity solver, i.e.
Hubbard-I, and the self energy $\hat{\Sigma}(i\omega)$ is then obtained, with
which we can construct the lattice Green’s function by applying the same self
energy term as,

$[\hat{G}^{\mathbf{k}}_{lattice}]^{-1}=\rmi\omega-\hat{H}_{\mathbf{k}}-\hat{\Sigma}_{\mathbf{k}}(\rmi\omega)+\mu$
(19)

The hybridization function, which characterizes the dynamics of the “heat
bath” is defined as
$\Delta(\rmi\omega)_{\alpha\beta}=\delta_{\alpha\beta}\sum\limits_{q}\frac{|V_{q\alpha}|^{2}}{\rmi\omega-\epsilon_{q\alpha}}$
and can be obtained iteratively by the following DMFT self consistent
equation.

$\hat{\Delta}(\rmi\omega)=\rmi\omega-\hat{E}_{imp}-\hat{\Sigma}+\mu-\left[\sum_{k}G^{\mathbf{k}}_{lattice}(\rmi\omega)\right]^{-1}$
(20)

The above equation is obtained by requiring that the local Green’s function on
the lattice should equal to the Green’s function of the quantum impurity
problem. The equations (19) and (20) form a closed self consistent loop, which
determines both the self energy and the hybridization function iteratively.

####  3.2.3 Double Counting Term

In this section we discuss the explicit expression of $E^{imp}_{\alpha}$ as
well as the double-counting term in (12). For 3d system under cubic symmetry,
the local interaction can be simplified as coulomb interaction $U$ and Hund’s
rule coupling $J$, which can be written as,

$\hat{H}_{U}=U\sum_{b}\hat{n}_{b,\uparrow}\hat{n}_{b,\downarrow}+(U-2J)\sum_{b<b^{\prime}\atop\sigma\sigma}\hat{n}_{b,\sigma}\hat{n}_{b^{\prime},\sigma}-J\sum_{b<b^{\prime}\atop\sigma}\hat{n}_{b\sigma}\hat{n}_{b^{\prime}\sigma}$
(21)

And the double counting energy in this case has been studied by Held et
al[21]. and can be approximately chosen as

$\displaystyle E_{DC}=\frac{1}{2}\bar{U}n_{\mbox{\tiny
LDA+DMFT}}(n_{\mbox{\tiny LDA+DMFT}}-1)$ (22)

where

$\displaystyle\bar{U}$ $\displaystyle=$
$\displaystyle\frac{U+2(M-1)(U-2J)-(M-1)J}{2M-1}$ (23) $\displaystyle M$
$\displaystyle=$ $\displaystyle 2l+1$ (24)

and $n_{\mbox{\tiny LDA+DMFT}}$ is the total number of electrons on correlated
orbitals for a specific atom, hence in (10) becomes

$\displaystyle\hat{V}_{DC}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\frac{\partial E_{DC}}{\partial
n_{\alpha}}\hat{n}_{\alpha}$ (25) $\displaystyle=$
$\displaystyle\sum_{\alpha}\bar{U}(n_{\mbox{\tiny
LDA+DMFT}}-\frac{1}{2})\hat{n}_{\alpha}$

The final expression of $E^{imp}$ is

$\displaystyle E^{imp}_{\alpha}$ $\displaystyle=$
$\displaystyle\sum_{k}\langle\alpha_{k}|\hat{H}_{LDA}|\alpha_{k}\rangle-
V_{DC}^{\alpha}$ (26)

where

$\displaystyle V_{DC}^{\alpha}=\bar{U}(n_{\mbox{\tiny LDA+DMFT}}-\frac{1}{2})$
(27)

The way to remove the form of double-counting energy is not unique, and in
fact this process needs intuition of physics. The double-counting discussed in
earlier sections is a very usual way which is also used in LDA+$U$ method.
Besides, the following two ways are also commonly used

  * •

to remove self-energy at zero-frequency
$\hat{\Sigma}(\rmi\omega)\rightarrow\hat{\Sigma}(\rmi\omega)-\hat{\Sigma}(0)$

  * •

to remove self-energy at infinity-frequency
$\hat{\Sigma}(\rmi\omega)\rightarrow\hat{\Sigma}(\rmi\omega)-\hat{\Sigma}(\infty)$

Furthermore, the double-counting energy is also regarded as “impurity solver
dependent”. For example, it is reasonable to use an integer to replace
$n_{\mbox{\tiny LDA+DMFT}}$ in (22) for the Hubbard-I solver[6], and in this
paper we also follow their suggestion.

####  3.2.4 Hubbard-I solver

The essential point of Hubbard-I solver is to neglect the effect of the heat
bath and take the atomic self energy as the zeroth order approximation for the
quantum impurity problem, which can be written as

$\displaystyle\hat{\Sigma}^{atom}=[\rmi\omega_{n}+\mu_{atom}-\hat{E}^{imp}]^{-1}-[\hat{G}^{atom}]^{-1}$
(28)

In the above equation $\hat{G}^{atom}$ is the Green’s function for a single
atom, which can be expressed as

$\displaystyle
G^{atom}_{\alpha\alpha^{\prime}}(\rmi\omega)=\sum_{\Gamma\Gamma^{\prime}}\frac{(F^{\alpha})_{\Gamma\Gamma^{\prime}}(F^{\alpha^{\prime}\dagger})_{\Gamma^{\prime}\Gamma}(X_{\Gamma}+X_{\Gamma^{\prime}})}{\rmi\omega-
E_{\Gamma^{\prime}}+E_{\Gamma}}$ (29)

where $|\Gamma\rangle$ and $|\Gamma^{\prime}\rangle$ are the atomic
eigenstates obtained by exact diagonalization of a single atom problem,
$X_{\Gamma^{\prime}}={e^{-\beta E_{\Gamma^{\prime}}}}/({\sum_{\Gamma}e^{-\beta
E_{\Gamma}}})$ represents the occupation probability of the local
configuration $|\Gamma\rangle$, and
$F^{\alpha}_{\Gamma\Gamma^{\prime}}=\langle\Gamma|f_{\alpha}|\Gamma^{\prime}\rangle$.

Hubbard-I approximation is

$\displaystyle\hat{\Sigma}\approx\hat{\Sigma}^{atom}$ (30)

In reference [22] , Savrasov et alpropose that the above self energy can be
written by the summation of a set of poles, which greatly simplifies the DMFT
calculation, because it is not necessary to handle the full frequency
dependence of the Green’s function. In Hubbard-I approximation, it can be
proved that both the atomic Green’s function and self energy are diagonal and
can be written in the form of pole expansion if the following two necessary
conditions are satisfied. 1) The single particle basis used here should be the
one which diagonalize the local density matrix; 2) The atomic Hamiltonian only
contains the two-body interaction terms. Obviously by chosen the proper single
particle basis the atomic Hamiltonian (18) in general will satisfy the above
two condition and thus can be always written in terms of pole expansion.
Unlike the reference[22], where they only use a few poles to capture the main
features of the self energy, in the present implementation, we keep all the
poles in the self energy, which makes it more accurate. From (29), it is very
clear that the atomic green’s function has already been expressed in terms of
poles and the pole expansion of the corresponding self energy can be obtained
using (28), which is introduced in details in the Appendix. Therefore in
general the above atomic self energy can be written as

$\hat{\Sigma}(\rmi\omega)=\hat{\Sigma}({\infty})+\sum^{np}_{i=1}\hat{V}^{\dagger}_{i}(\rmi\omega-
p_{i})^{-1}\hat{V}_{i}$ (31)

where $i$ labels the number of poles, $np$ is the total number of poles in the
self energy and $\hat{V}_{i}$ is a vector defined in the orbital space
describing the distribution of the $ith$ “pole states” among the local
orbitals.

###  3.3 Correction of the Density Matrices and Pole Expansion of the Self-
energy

Once we obtain the converged local self-energy, we are at the point to correct
the k-dependent lattice Green’s functions as well as the density matrices. In
general we need to do the summation over the Matsubara frequencies, which is
time consuming for realistic materials contains lots of bands. But it can be
greatly simplified when the self energy can be expressed as the summation of
poles, which can be written as

$\displaystyle\langle\varphi_{n_{1}\mathbf{k}}|\hat{\rho}_{LDA+DMFT}|\varphi_{n_{2}\mathbf{k}}\rangle$
(32) $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{i\omega_{n}}\langle\varphi_{n_{1}\mathbf{k}}|\frac{1}{\rmi\omega_{n}+\mu-\hat{H}_{LDA}+\hat{V}_{DC}-\hat{\Sigma}(i\omega_{n})}|\varphi_{n_{2}\mathbf{k}}\rangle$

As pointed out firstly in reference [22], when the self energy can be written
in terms of poles, the full green’s function can be expressed as the
“physical” part of an enlarged “pseudo Hamiltonian”, defined as

$\displaystyle\hat{H}_{ps}(\mathbf{k})=\left(\begin{array}[]{cc}\hat{H}_{LDA}(\mathbf{k})-\hat{V}_{DC}+\hat{\Sigma}(\infty)&\qquad-\hat{V}^{\dagger}\\\
-\hat{V}&\qquad\hat{P}\end{array}\right)$ (35)

where
$\hat{H}_{LDA}(\mathbf{k})=E_{n\mathbf{k}}\hat{c}^{\dagger}_{n\mathbf{k}}\hat{c}_{n\mathbf{k}}$
is the physical Hilbert space, $\hat{V}_{i}$ is defined in (31) and
$\hat{P}=Diag(p_{1},p_{2},....,p_{np})$. Therefore the physical green’s
function can be expressed in terms of eigenstate of the pseudo Hamiltonian
simply as

$\displaystyle\langle\varphi_{n_{1}\mathbf{k}}\bigg{|}\frac{1}{\rmi\omega_{n}-\hat{H}_{LDA}(\mathbf{k})+\hat{V}_{DC}-\hat{\Sigma}(i\omega_{n})}\bigg{|}\varphi_{n_{2}\mathbf{k}}\rangle$
(36) $\displaystyle=$
$\displaystyle\langle\varphi_{n_{1}\mathbf{k}}|\psi^{ps}_{l\mathbf{k}}\rangle\frac{1}{\rmi\omega_{n}-E^{ps}_{l}(\mathbf{k})}\langle\psi^{ps}_{l\mathbf{k}}|\varphi_{n_{2}\mathbf{k}}\rangle$

where $|\psi^{ps}_{l\mathbf{k}}\rangle$ and $E^{ps}_{l}(\mathbf{k})$ are the
eigenstate and eigenvalue of pseudo Hamiltonian respectively.

Then the sum of frequencies in (32) can be performed directly

$\displaystyle\langle\varphi_{n_{1}\mathbf{k}}|\hat{\rho}_{LDA+DMFT}|\varphi_{n_{2}\mathbf{k}}\rangle$
(37) $\displaystyle=$
$\displaystyle\frac{1}{\beta}\sum_{i\omega_{n}}\sum_{l\mathbf{k}}\langle\varphi_{n_{1}\mathbf{k}}|\psi^{ps}_{l\mathbf{k}}\rangle\frac{1}{\rmi\omega_{n}-E_{l}^{ps}(\mathbf{k})}\langle\psi^{ps}_{l\mathbf{k}}|\varphi_{n_{2}\mathbf{k}}\rangle$
$\displaystyle=$
$\displaystyle\sum_{l\mathbf{k}}\langle\varphi_{n_{1}\mathbf{k}}|\psi^{ps}_{l\mathbf{k}}\rangle\langle\psi^{ps}_{l\mathbf{k}}|\varphi_{n_{2}\mathbf{k}}\rangle
n_{F}[E_{l}^{ps}(\mathbf{k})-\mu]$

in which $\mu$ is exactly the chemical potential in (12), determined by the
occupation number of electrons $N_{tot}$ in the unit cell

$\displaystyle N_{tot}$ $\displaystyle=$
$\displaystyle\sum_{n\mathbf{k}}\langle\varphi_{n\mathbf{k}}|\hat{\rho}|\varphi_{n\mathbf{k}}\rangle$
(38) $\displaystyle=$
$\displaystyle\sum_{ln\mathbf{k}}|\langle\varphi_{n\mathbf{k}}|\psi^{ps}_{l\mathbf{k}}\rangle|^{2}n_{F}[E_{l}^{ps}(\mathbf{k})-\mu]$

###  3.4 Evaluation of Real Space Charge Density and Complete of the Full Loop

After DMFT obtains the corrected density matrix $\hat{\rho}$, the real space
charge density will be generated based on it by Fourier transformation, and
then, the new LDA Hamiltonian as well as the new overlap matrices between
Bloch states and local orbitals will be recalculated, which closes the full
iteration loop for the charge density self consistency. The occupation number
of the Bloch basis can also be obtained by,

$\displaystyle n_{\mbox{\tiny
LDA+DMFT}}=\sum_{\mathbf{k},n_{1}n_{2},\alpha}\langle\alpha_{\mathbf{k}}|\varphi_{n_{1}\mathbf{k}}\rangle\langle\varphi_{n_{1}\mathbf{k}}|\hat{\rho}|\varphi_{n_{2}\mathbf{k}}\rangle\langle\varphi_{n_{2}\mathbf{k}}|\alpha_{\mathbf{k}}\rangle$
(39)

###  3.5 Calculation of Physical Quantities

####  3.5.1 Energy Functional

Most of the physical quantities we are interested in can be calculated once we
have reached the charge density self consistency. First we present the energy
functional of the full loop LDA+DMFT, which can be used to calculate many
quantities, i.e. the total energy, force and so on. It can be written as, [6]

$\displaystyle E=E_{LDA}[\rho]-\langle H_{KS}\rangle_{LDA}+\langle
H_{KS}\rangle+\langle H_{U}\rangle-E_{DC}$ (40)

in which $E_{LDA}[\rho]$ is the expression of the energy within density-
functional theory; $\langle
H_{KS}\rangle_{LDA}=\tr[\hat{H}_{LDA}\hat{G}_{LDA}]$ is the non-interacting
energy at LDA level; $\langle H_{KS}\rangle=\tr[\hat{H}_{LDA}\hat{G}]$ is the
non-interacting energy at DMFT level; $\langle
H_{U}\rangle=\frac{1}{2}\tr[\hat{\Sigma}\hat{G}]$ is the interaction energy
caused by local correlation interactions; the last term $E_{DC}$ is the double
counting energy given before in (22). In above the meaning for “trace” is
defined as

$\displaystyle\tr[\mathcal{A}]=\frac{1}{\beta}\sum_{n\mathbf{k}}\sum_{\rmi\omega_{n}}\langle\varphi_{n\mathbf{k}}|\mathcal{A}(\rmi\omega_{n})|\varphi_{n\mathbf{k}}\rangle$
(41)

and the integral path surrounds all the energies of occupation states.

The first and the second term is easily calculated in the LDA framework, and
the three terms remaining must be evaluated in the DMFT process, explicitly

$\displaystyle\langle
H_{KS}\rangle=\sum_{n\mathbf{k}}E_{n\mathbf{k}}\langle\varphi_{n\mathbf{k}}|\hat{\rho}|\varphi_{n\mathbf{k}}\rangle$
(42)

and

$\displaystyle\langle H_{U}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\tr[\hat{\Sigma}\hat{G}]$ (43) $\displaystyle=$
$\displaystyle\frac{1}{2}\tr[(\hat{G}^{-1}_{0}-\hat{G}^{-1})\hat{G}]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\sum_{m,n\mathbf{k}}(E^{ps}_{m\mathbf{k}}-\mu)|\langle\varphi_{n\mathbf{k}}|\psi_{m\mathbf{k}}\rangle|^{2}n_{F}(E^{ps}_{m\mathbf{k}}-\mu)\right.$
$\displaystyle\left.-\langle H_{KS}\rangle+\langle\hat{V}_{DC}\rangle\right]$

while $\langle H_{KS}\rangle$ is expressed in (42), and
$\langle\hat{V}_{DC}\rangle$ is the expected value of $\hat{V}_{DC}$ at DMFT
level

$\displaystyle\langle\hat{V}_{DC}\rangle$ $\displaystyle=$
$\displaystyle\sum_{i}\bar{U}_{i}(n^{i}_{\mbox{\tiny
LDA+DMFT}}-\frac{1}{2})\sum_{\alpha}\langle\hat{n}^{i}_{\alpha}\rangle$ (44)
$\displaystyle=$ $\displaystyle\bar{U}_{i}(n^{i}_{\mbox{\tiny
LDA+DMFT}}-\frac{1}{2})n^{i}_{\mbox{\tiny LDA+DMFT}}$ (45)

Notice that here emerges the site index $i$. We should be very careful in the
case that the number of correlated atoms is larger than one. The final
expression of the total energy is

$\displaystyle E=E_{LDA}[\rho]-\langle H_{KS}\rangle_{LDA}+\langle
H_{KS}\rangle$
$\displaystyle+\frac{1}{2}\left[\sum_{m,nk}(E^{ps}_{mk}-\mu)|\langle\varphi_{nk}|\psi_{mk}\rangle|^{2}n_{F}(E^{ps}_{mk}-\mu)-\langle
H_{KS}\rangle\right]$
$\displaystyle+\frac{1}{4}\sum_{i}\bar{U}_{i}n^{i}_{\mbox{\tiny LDA+DMFT}}$
(46)

However, if we use Hubbard-I as the impurity solver, and remove double-
counting term by using an integer, the energy functional is given by Haule[7]

$\displaystyle E=E_{LDA}[\rho]-\langle H_{KS}\rangle_{LDA}+\langle
H_{KS}\rangle$
$\displaystyle+\frac{1}{2}\left[\sum_{m,n\mathbf{k}}(E^{ps}_{m\mathbf{k}}-\mu)|\langle\varphi_{n\mathbf{k}}|\psi_{m\mathbf{k}}\rangle|^{2}n_{F}(E^{ps}_{m\mathbf{k}}-\mu)-\langle
H_{KS}\rangle\right]$ $\displaystyle+\frac{1}{4}\sum_{i}\bar{U}_{i}n^{i}_{0}$
(47)

where $n^{i}_{0}$ refers to the integer used to remove double counting for the
$i$th correlated atom.

####  3.5.2 Expressions of DOS and PDOS

In order to calculate the density of state (DOS) or the partial density of
state (PDOS) on correlated orbitals, it is necessary to calculate the retarded
form of the Green’s function. By using the virtue of pole expansion method, we
write directly

$\displaystyle
D(\epsilon)=(-\frac{1}{\pi})\mathrm{Im}\bigg{[}\sum_{\mathbf{k},m,n}\frac{\langle\varphi_{n\mathbf{k}}|\psi_{m\mathbf{k}}\rangle\langle\psi_{m\mathbf{k}}|\varphi_{n\mathbf{k}}\rangle}{\epsilon+\mu-E^{ps}_{m\mathbf{k}}+\rmi\eta}\bigg{]}\
(\eta\rightarrow 0^{+})$ (48)

and the PDOS of orbital $\alpha$

$\displaystyle D_{\alpha}(\epsilon)$ $\displaystyle=$
$\displaystyle(-\frac{1}{\pi})\mathrm{Im}\left[\sum_{\mathbf{k},m}\frac{\langle\alpha_{\mathbf{k}}|\psi_{m\mathbf{k}}\rangle\langle\psi_{m\mathbf{k}}|\alpha_{\mathbf{k}}\rangle}{\epsilon+\mu-E^{ps}_{m\mathbf{k}}+\rmi\eta}\right]\
(\eta\rightarrow 0^{+})$ (49)

where

$\displaystyle\langle\alpha_{\mathbf{k}}|\psi_{m\mathbf{k}}\rangle=\sum_{n}\langle\alpha_{\mathbf{k}}|\varphi_{n\mathbf{k}}\rangle\langle\varphi_{n\mathbf{k}}|\psi_{m\mathbf{k}}\rangle$
(50)

##  4 Benchmark

To validate the LDA+DMFT framework based on pseudo-potentials-planewave
package and pole expansion of self energy, we benchmark our implementation by
applying to $\gamma$-cerium, americium and paramagnetic NiO. These three are
typical strongly correlated systems in which the valence electrons are
believed to be on the localized side as reported in previous literatures. We
use the Hubbard-I method as impurity solver, which is good enough to capture
the atomic-like features in the Mott insulators, so we expect that these
systems could be well described by our method. We would emphasize here that in
our implementation we can replace the Hubbard-I solver by any of the solvers
as long as the self energy can be written in pole expansion form.

###  4.1 Cerium

The cerium metal attracts lots of research interests for its isostructural
volume-collapse transition from $\gamma$ phase to $\alpha$ phase. The volume
change is about $15\%$ during the transition, which is possibly driven by the
entropy change[23]. Normal LDA calculations could not give a correct
description of cerium, as show in table 1, the equillibrium volume of cerium
given by LDA calculations is smaller than the experimental volume of $\alpha$
phase. This is due to the fact that in LDA the $f$-electrons are treated as
itinerant and the strong correlation effects among them can not be well
captured. When $f$-electron is treated as core electrons, we see that the
equilibrium volume is larger but close to $\gamma$ phase. These simple LDA
results implies that in $\gamma$ phase, the $f$-electrons are more localized
than in $\alpha$-phase. Therefore $f$-electrons in $\gamma$ phase is quite
close to the localized picture like the situation in Mott insulators, which
makes it suitable to applying the Hubbard-I solver. The $\alpha$-Ce which is
stable in low temperature is a correlated metal, which is confirmed by many
other LDA+DMFT calculations[23], and also by our newly development
LDA+Gutzwiller method[24].

In our calculation for $\gamma$-Ce, we used ultrasoft pseudo-potentional in
which 4 _f_ , 5 _s_ , 5 _d_ , 6 _s_ orbitals are taken as valence orbitals.
The energy cutoff for plane-wave expansion is 12.5 Ha for convergence.
Calculations are performed with $10\times 10\times 10$ k-points. For DMFT
calculation, the on-site Coulomb interaction $U$ is chosen to be 6.0 eV as
suggested in previous reports.[25, 23, 26]

In figure 2 the total energy versus volume for $\gamma$-cerium obtained by
LDA, LDA+$U$ and LDA+DMFT are plotted. The equilibrium volume and bulk modulus
are also shown in table 1. We can find in figure 2 that the equilibrium volume
obtained by our LDA+DMFT calculation is very close to the experimental data,
which shows great improvements over LDA . The bulk modulus is also improved
but still larger than the experimental data, which may imply the contribution
for lattice vibration is inneglectable as discussed in reference [6].

Table 1:  Lattice parameter and bulk modulus of $\gamma$-cerium according to
both experiment and calculations

  * | Volume($\mathrm{\AA}^{3}$) | Bulk modulus(GPa)
---|---|---
Experiments | 34.35[27] | 21[28]
LDA | 22.36 | 62.09
LDA(fcore)[29] | 36.39 | 30.2
LDA+$U$/PAW[30] | 32.0 | 34
LDA+DMFT[6] | 30.10 | 48.49
LDA+DMFT | 33.29 | 38.27

Figure 2: Total energy versus volume curves obtained by LDA, LDA+$U$[6], and
LDA+DMFT. The dashed lines are fitted by Birch EOS. The energy were shifted
for different curves. The experimental equilibrium volume is from reference
[27]

Figure 3 shows the 4 _f_ partial DOS for $\gamma$ cerium at the equilibrium
volume, the position of the two Hubbard bands agree well with the experimental
result.

Figure 3: 4 _f_ pdos for $\gamma$ cerium in LDA+DMFT (Hubbard I), XPS and BIS
experimental data are from reference [31] and [32]

###  4.2 Americium

Americium, which is widely used in smoke detectors, can be viewed as the
”changing point” of the actinide series, where the 5 _f_ electrons changes
from delocalized to localized states[33]. Valence-band ultraviolet
photoemission experiment by Naegele et al[34] shows that the 5 _f_ states are
strongly localized. Am undergoes a series of structure phase transitions with
the increment of pressure, from double hexagonal close packed(dhcp, P63/mmc),
to face centered cubic(fcc, Fm3m, 6.1 GPa), face centered orthorhombic(Fddd,
10.0 GPa), and to primitive orthorhombic structure(Pnma, 16 PGa)[35].
Therefore Am provides an interesting playfield to investigate the relation
between $f$-electron localization and structure transition. Here we focus on
the fcc phase and dhcp phase.

In the calculation, we used ultrasoft pseudo-potentional contains 5 _f_ , 6
_p_ , 6 _d_ , 7 _s_ orbitals as valence orbitals, the plane-wave expansion
kinetic energy cutoff of 12.5 Ha. The LDA self-consistent calculations were
performed with $10\times 10\times 10$ k-points grid for fcc phase, $9\times
9\times 3$ k-points grid for dhcp phase. For the volume calculation, we only
considered density-density interaction, and on-site Coulomb interaction
$F^{(0)}=4.5$ eV, for the partial DOS calculation, we considered general
interaction, and take the atomic values $F^{(2)}=7.2$ eV, $F^{(4)}=4.8$ eV,
$F^{(6)}=3.6$ eV[36].

We present the optimized volume and bulk modulus for both dhcp and fcc phase
in table 2 and figure 4. Our LDA+DMFT results are quite close to the
experimental data and show great improvement over LDA. Also our results are
quite consistent with the reference [22], in which the full potential LMTO
methods are used in the LDA part. The equation of state of Am obtained by our
LDA+DMFT calculation is plotted in figure 5, which again shows very good
agreement with both the experimental data and previous LDA+DMFT calculation
based on LMTO method[22].

Table 2: Am equilibrium volume and bulk modulus for dhcp and fcc phase by
different methods, together with the experimental results[35, 37].

  * | Volume($\mathrm{\AA}^{3}$) | Bulk modulus(GPa)
---|---|---
experiment | 29.250 | 29.9
GGA(dhcp)[38] | 19.916 | 70.0
LDA(dhcp) | 18.55 | 118.99
LDA(fcc) | 17.73 | 169.40
LDA+DMFT(dhcp) | 27.04 | 52.70
LDA+DMFT(fcc) | 26.99 | 50.75
LDA+DMFT[22] | 27.4 | 45

Figure 4: Total energy per atom of Am obtained by LDA and LDA+DMFT. The
LDA+DMFT results of Am for both phase are very close, and the experimental
equilibrium volume is denoted by dot line. Figure 5: Calculated P-V relation
of dhcp and fcc Am. Here, $V_{0}$ is experimental equilibrium volume, The
experimental P-V results[37] are shown by black dots.

Figure 6 shows the 5 _f_ partial DOS of fcc phase Am obtained both by our
LDA+DMFT calculation and LDA+DMFT(OCA) by Savrasov[22],together with the
photoemission data from ref.[34]. Our results is quite consistent with the
previous calculations and the experimental results.

Figure 6: The spectral function of Am evaluated by LDA+DMFT. The experimental
photoemission results is from reference [34] are denoted by solid dots.

###  4.3 Paramagnetic phase of NiO

NiO has been heavily studied as a prototype of Mott insulators. It has been
studied within the frame of LDA+DMFT by many groups[39, 40]. Therefore NiO can
be used as a good benchmark material for the implementation of LDA+DMFT. In
this paper, we apply our code to study the paramagnetic phase of NiO. First we
plot the partial density of states (PDOS) obtained by LDA in figure 7, which
clearly shows metallic behavior and is not consistent with the experiments. In
our LDA+DMFT(Hub1) calculation, we use the wannier functions as the local
basis and choose on-site interaction $U=8.0eV$, Hund’s coupling $J=1.0eV$. The
Mott insulating features of NiO then can be well captured by our LDA+DMFT
method, as shown by the PDOS in figure 8).

Figure 7: The partial DOS of NiO obtained by LDA. The main contribution to the
Fermi surface is attributed to eg-like Wannier orbitals, and the t2g-like
orbitals are almost fully occupied. Figure 8: The spectral function of NiO
obtained by LDA+DMFT calculation. As shown in the figure, compared with LDA
results, the t2g-like orbitals are still fully occupied while the eg-like
orbitals split into upper and lower Hubbard bands.

The energy gap obtained by our calculation is around 4.0eV, which is also
quite consistent with the experiments. The photo emission and BIS data
obtained by Sawatzky and Allen[41] are also plotted in figure 9. We can find
that our results fit the photo emission and BIS data very well.

Figure 9: The spectral function of NiO evaluated by LDA+DMFT, with the
comparison of experimental data[41]. The behavior of the density of state near
the fermi surface fits well with the experimental data.

##  5 Conclusions

The new implementation of LDA+DMFT based on the pseudo-potential plane-wave
method is introduced in detail in this paper. We choose the Hubbard-I method
as the impurity solver to solve the quantum impurity problem generated by
DMFT, which is quite suitable for the Mott insulator materials. The most
important advantage of Hubbard-I method is that it is simple enough, which
makes the full loop charge self consistent calculation accessible. We also
point out in the paper that the difficulty of handling frequency dependent
Green’s function can be completely removed by expressing the self energy in
terms of pole expansion, which greatly raise the efficiency of the method.
Finally we benchmark our implementation of LDA+DMFT on several important
correlation materials including $\gamma$-Ce,Am and NiO. Our results for all
these materials fit very well with both the experimental data and the previous
LDA+DMFT results.

##  Appendix A Self-energy in pole expansion form

In this appendix, we will introduce how to evaluate self-energy from Green’s
function in the pole expansion form.

Green’s function can be expressed as

$G(i\omega)=\sum^{N_{G}}_{i=1}\frac{V_{i}}{i\omega-P_{i}}$ (51)

where $P_{i}$ is the $i$th pole of the Green’s function, $V_{i}$ is the weight
of the pole, and $N_{G}$ is the number of Green’s function poles. Besides the
self-energy can be expressed as

$\Sigma(i\omega)=\sum^{N_{S}}_{i=1}\frac{W_{i}}{i\omega-Q_{i}}+\Sigma(\infty)$
(52)

where $Q_{i}$ is the $i$th pole of the self-energy, $W_{i}$ is the weight of
the pole, and $N_{S}$ is the number of self-energy poles.

Then the Dyson’s equation (28) becomes

$\sum^{N_{S}}_{i=1}\frac{W_{i}}{i\omega-
Q_{i}}+\Sigma(\infty)=i\omega+\mu-\epsilon_{imp}-\left[\sum^{N_{G}}_{i=1}\frac{V_{i}}{i\omega-
P_{i}}\right]^{-1}$ (53)

First consider the $\Sigma(\infty)$, when $\omega\rightarrow\infty$

$\displaystyle\Sigma(\omega\rightarrow\infty)$
$\displaystyle=i\omega+\mu-\epsilon_{imp}-\left[\sum^{N_{G}}_{i=1}\frac{V_{i}}{i\omega-
P_{i}}\right]^{-1}$ (54)
$\displaystyle=i\omega+\mu-\epsilon_{imp}-\left[\frac{1}{i\omega}\sum^{N_{G}}_{i=1}\frac{V_{i}}{1-\frac{P_{i}}{i\omega}}\right]^{-1}$
$\displaystyle=i\omega+\mu-\epsilon_{imp}-i\omega\left[\sum^{N_{G}}_{i=1}V_{i}-\frac{1}{i\omega}\sum^{N_{G}}_{i=1}V_{i}P_{i}\right]^{-1}$
$\displaystyle=i\omega+\mu-\epsilon_{imp}-i\omega\left(1-\frac{\sum^{N_{G}}_{i=1}V_{i}P_{i}}{i\omega}\right)$
$\displaystyle=\mu-\epsilon_{imp}+\sum^{N_{G}}_{i=1}V_{i}P_{i}$

For the the self energy poles in (53), it corresponding to the zero in the
square bracket of the right side of (53), that is for $Q\in\\{Q_{i}\\}$, we
have

$\sum^{N_{G}}_{i=1}\frac{V_{i}}{Q-P_{i}}=0$ (55)

The left hand side is monotonically decreasing between two adjacent Green’s
function poles, so we could use bisection method to find the self energy
poles.

We rewrite (53) as

$\sum^{N_{s}}_{i=1}\frac{W_{i}}{i\omega-
Q_{i}}=i\omega+\mu-\epsilon-\Sigma(\infty)-\left[\sum^{N_{g}}_{i=1}\frac{V_{i}}{i\omega-
P_{i}}\right]^{-1}$ (56)

For the left hand side

$\displaystyle W_{i}$ $\displaystyle=\lim_{i\omega\rightarrow
Q_{i}}\left(\sum^{N_{S}}_{j=1}\frac{W_{j}}{i\omega-Q_{j}}\right)(i\omega-
Q_{i})$ (57) $\displaystyle=\lim_{i\omega\rightarrow
Q_{i}}\left[i\omega+\mu-\epsilon_{imp}-\Sigma(\infty)-\left(\sum^{N_{g}}_{j=1}\frac{V_{j}}{i\omega-
P_{j}}\right)^{-1}\right](i\omega-Q_{i})$
$\displaystyle=\lim_{i\omega\rightarrow Q_{i}}-(i\omega-
Q_{i})\left(\sum^{N_{G}}_{j=1}\frac{V_{j}}{i\omega-P_{j}}\right)^{-1}$

here we define $\delta=i\omega-Q_{i}$ and
$h(\delta)=\sum^{N_{G}}_{j=1}\frac{V_{j}}{Q_{i}+\delta-P_{j}}$, expand
$h(\delta)$ as

$h(\delta)\approx\delta\cdot
h^{\prime}(0)=-\sum^{N_{G}}_{j=1}\frac{V_{j}}{(Q_{i}-P_{j})^{2}}\cdot\delta$
(58)

so

$W_{i}=\left[\sum^{N_{G}}_{j=1}\frac{V_{j}}{(Q_{i}-P_{j})^{2}}\right]^{-1}$
(59)

here we have got both self energy poles and its weight.

The authors thank L.Wang, G. Xu, and F. Lu for their helpful discussions. We
acknowledge the support from the 973 Program of China (Grant No.
2011CBA00108), from the NSF of China (Grants No. NSFC 10876042 and No. NSFC
10874158), and that from the Fund of Laboratory of Shock Wave and Detonation
Physics, Institute of Fluid Physics, CAEP, Contract No. 2011-056000-0833F.

## References

## References

  * [1] Georges A, Kotliar G, Krauth W, and Rozenberg M J.  Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions.  Rev. Mod. Phys., 68(1):13–125, January 1996.
  * [2] Lechermann F, Georges A, Poteryaev A, Biermann S, Posternak M, Yamasaki A, and Andersen O K.  Dynamical mean-field theory using Wannier functions: A flexible route to electronic structure calculations of strongly correlated materials.  Phys. Rev.B, 74(12):125120, September 2006.
  * [3] Anisimov V I, Kozhevnikov A V, Korotin M A, Lukoyanov A V, and Khafizullin D A.  Orbital density functional as a means to restore the discontinuities in the total-energy derivative and the exchange–correlation potential.  J. Phys.: Condens. Matter, 19(10):106206, March 2007.
  * [4] Savrasov S Y and Kotliar G.  Spectral density functionals for electronic structure calculations.  Phys. Rev.B, 69(24):245101, June 2004.
  * [5] Minár J, Chioncel L, Perlov A, Ebert H, Katsnelson M, and Lichtenstein A.  Multiple-scattering formalism for correlated systems: A KKR-DMFT approach.  Phys. Rev.B, 72(4):045125, July 2005.
  * [6] Pourovskii L V, Amadon B, Biermann S, and Georges A.  Self-consistency over the charge density in dynamical mean-field theory: A linear muffin-tin implementation and some physical implications.  Phys. Rev.B, 76(23):235101, December 2007.
  * [7] Haule K, Yee Chuck-Hou, and Kim Kyoo.  Dynamical mean-field theory within the full-potential methods: Electronic structure of CeIrIn5 , CeCoIn5 , and CeRhIn5.  Phys. Rev.B, 81(19):195107, May 2010.
  * [8] Andersen O K.  Linear methods in band theory.  Phys. Rev.B, 12(8):3060–3083, October 1975.
  * [9] Anisimov V I, Aryasetiawan F, and Lichtenstein A I.  First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA + U method.  J. Phys.: Condens. Matter, 9(4):767–808, January 1997.
  * [10] Lichtenstein A I and Katsnelson M I.  Ab initio calculations of quasiparticle band structure in correlated systems: LDA++ approach.  Phys. Rev.B, 57(12):6884–6895, March 1998.
  * [11] Savrasov S Y, Kotliar G, and Abrahams E.  Correlated electrons in delta-plutonium within a dynamical mean-field picture.  Nature, 410(6830):793–5, April 2001.
  * [12] Andersen O K and Saha-Dasgupta T.  Muffin-tin orbitals of arbitrary order.  Phys. Rev.B, 62(24):R16219–R16222, December 2000.
  * [13] Anisimov V I, Kondakov D E, Kozhevnikov A V, Nekrasov I A, Pchelkina Z V, Allen J W, Mo S K, Kim H D, Metcalf P, Suga S, Sekiyama A, Keller G, Leonov I, Ren X, and Vollhardt D.  Full orbital calculation scheme for materials with strongly correlated electrons.  Phys. Rev.B, 71(12):125119, March 2005.
  * [14] Souza I, Marzari N, and Vanderbilt D.  Maximally localized Wannier functions for entangled energy bands.  Phys. Rev.B, 65(3):035109, December 2001.
  * [15] Marzari N and Vanderbilt D.  Maximally localized generalized Wannier functions for composite energy bands.  Phys. Rev.B, 56(20):12847–12865, November 1997.
  * [16] Amadon B, Lechermann F, Georges A, Jollet F, Wehling T O, and Lichtenstein A I.  Plane-wave based electronic structure calculations for correlated materials using dynamical mean-field theory and projected local orbitals.  Phys. Rev.B, 77(20):205112, May 2008.
  * [17] Trimarchi G, Leonov I, Binggeli N, Korotin Dm, and Anisimov V I.  LDA+DMFT implemented with the pseudopotential plane-wave approach.  J. Phys.: Condens. Matter, 20(13):135227, April 2008.
  * [18] Korotin DM, Kozhevnikov A V, Skornyakov S L, Leonov I, Binggeli N, Anisimov V I, and Trimarchi G.  Construction and solution of a Wannier-functions based Hamiltonian in the pseudopotential plane-wave framework for strongly correlated materials.  Eur. Phys. J. B, 98:91–98, 2008.
  * [19] Vanderbilt D.  Soft self-consistent pseudopotentials in a generalized eigenvalue formalism.  Phys. Rev.B, 41(11):7892–7895, April 1990.
  * [20] Hamann D R, Schlüter M, and Chiang C.  Norm-Conserving Pseudopotentials.  Phys. Rev. Lett., 43(20):1494–1497, November 1979.
  * [21] Held K.  Electronic structure calculations using dynamical mean field theory.  Advances in Physics, 56(6):829–926, November 2007.
  * [22] Savrasov S Y, Haule K, and Kotliar G.  Many-Body Electronic Structure of Americium Metal.  Phys. Rev. Lett., 96(3):036404, January 2006.
  * [23] Amadon B, Biermann S, Georges A, and Aryasetiawan F.  The $\alpha$-$\gamma$ Transition of Cerium Is Entropy Driven.  Phys. Rev. Lett., 96(6):066402, February 2006.
  * [24] Tian M F, Deng X Y, Fang Z, and Dai X.  The intermediate pressure phases of cerium: a LDA+Gutzwiller method study.  arXiv, page 1104.0156, April 2011.
  * [25] Haule K, Oudovenko V, Savrasov S Y, and Kotliar G.  The $\alpha\rightarrow\gamma$ transition in ce: A theoretical view from optical spectroscopy.  Phys. Rev. Lett., 94(3):036401, Jan 2005.
  * [26] Amadon B.  Fully self-consistent LDA+DMFT Calculation: A plane wave and projector augmented wave implementation and application to Cerium, Ce2O3 and Pu2O3.  arXiv, page 1101.0539, January 2011.
  * [27] Jeong I K, Darling T W, Graf M J, Proffen Th, Heffner R H, Lee Y, Vogt T, and Jorgensen J D.  Role of the lattice in the $\gamma\rightarrow\alpha$ phase transition of ce: A high-pressure neutron and x-ray diffraction study.  Phys. Rev. Lett., 92(10):105702, Mar 2004.
  * [28] Decremps F, Antonangeli D, Amadon B, and Schmerber G.  Role of the lattice in the two-step evolution of $\gamma$-cerium under pressure.  Phys. Rev.B, 80(13):132103, October 2009.
  * [29] Huang L and Chen C A.  The lattice dynamics of $\alpha$\- and $\gamma$-Ce: a first-principles approach.  J. Phys.: Condens. Matter, 19(47):476206, November 2007.
  * [30] Amadon B, Jollet F, and Torrent M.  $\gamma$ and $\beta$ cerium: LDA+U calculations of ground-state parameters.  Phys. Rev.B, 77(15):155104, April 2008.
  * [31] Wuilloud E, Moser H R, Schneider W D, and Baer Y.  Electronic structure of $\gamma$\- and $\alpha$-ce.  Phys. Rev.B, 28(12):7354–7357, Dec 1983.
  * [32] Wieliczka D, Weaver J H, Lynch D W, and Olson C G.  Photoemission studies of the $\gamma$-$\alpha$ phase transition in ce: Changes in 4$f$ character.  Phys. Rev.B, 26(12):7056–7059, Dec 1982.
  * [33] Moore K T and van der Laan G.  Nature of the 5f states in actinide metals.  Rev. Mod. Phys., 81(1):235–298, February 2009.
  * [34] Naegele J R, Manes L, Spirlet J C, and Müller W.  Localization of 5f Electrons in Americium: A Photoemission Study.  Phys. Rev. Lett., 52(20):1834–1837, May 1984.
  * [35] Heathman S, Haire R G, Le Bihan T, Lindbaum A, Litfin K, Méresse Y, and Libotte H.  Pressure Induces Major Changes in the Nature of Americium’s 5f Electrons.  Phys. Rev. Lett., 85(14):2961–2964, October 2000.
  * [36] Carnall W T and Wybourne B G.  Electronic Energy Levels of the Lighter Actinides: U3+, Np3+, Pu3+, Am3+, and Cm3+.  J. Chem. Phys., 40(11):3428, 1964.
  * [37] Lindbaum A, Heathman S, Litfin K, Méresse Y, Haire R G, Le Bihan T, and Libotte H.  High-pressure studies of americium metal: Insights into its position in the actinide series.  Phys. Rev.B, 63(21):214101, April 2001.
  * [38] Kutepov A and Kutepova S.  First-principles study of electronic and magnetic structure of alpha-Pu, delta-Pu, americium, and curium.  J. Magn. Magn. Mater., 272-276:E329–E330, May 2004.
  * [39] Ren X, Leonov I, Keller G, Kollar M, Nekrasov I, and Vollhardt D.  LDA+DMFT computation of the electronic spectrum of NiO.  Phys. Rev.B, 74(19):195114, November 2006.
  * [40] Kuneš J, Anisimov V I, Lukoyanov A V, and Vollhardt D.  Local correlations and hole doping in NiO: A dynamical mean-field study.  Phys. Rev.B, 75(16):4–7, April 2007.
  * [41] Sawatzky G and Allen J.  Magnitude and Origin of the Band Gap in NiO.  Phys. Rev. Lett., 53(24):2339–2342, December 1984.


# Towards Fast and Accurate Multi-Person Pose Estimation on Mobile Devices

Xuan Shen1 These Authors contributed equally. Geng Yuan∗1 Wei Niu2 Xiaolong
Ma1
Jiexiong Guan2 Zhengang Li1 Bin Ren2&Yanzhi Wang1 1Northeastern University,
2College of William & Mary
{shen.xu, yuan.geng, ma.xiaol, li.zhen<EMAIL_ADDRESS>
{wniu<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The rapid development of autonomous driving, abnormal behavior detection, and
behavior recognition makes an increasing demand for multi-person pose
estimation-based applications, especially on mobile platforms. However, to
achieve high accuracy, state-of-the-art methods tend to have a large model
size and complex post-processing algorithm, which costs intense computation
and long end-to-end latency. To solve this problem, we propose an architecture
optimization and weight pruning framework to accelerate inference of multi-
person pose estimation on mobile devices. With our optimization framework, we
achieve up to 2.51$\times$ faster model inference speed with higher accuracy
compared to representative lightweight multi-person pose estimator.

##  1 Introduction

In recent years, with the increasing popularity of mobile computing and edge
AI, mobile platforms have become an important carrier of DNN applications Li
et al. (2021a); Zhao et al. (2020); Niu et al. (2020). multi-person pose
estimation (PE) is one of the popular DNN applications that play an important
role in many fields, including autonomous driving, abnormal behavior
detection, and behavior recognition, and has made impressive progress with the
joint efforts.

The multi-person PE can be achieved mainly in two design methodologies, which
are the top-down method and the bottom-up method. The top-down method Fang et
al. (2017) performs PE on a single person who is captured by a detector in
advance. The accuracy will be dependent on the detector which might fail to
find all the persons when they are close together. And the overall inference
time will increase with the increasing number of people in an image frame,
which is not efficient for multi-person PE. On the other hand, the bottom-up
method Cao et al. (2019) predicts the human keypoints for all people in an
image frame at the same time. The bottom-up methods usually have smaller model
sizes and fewer computation counts (i.e., GFLOPs), making the bottom-up
methods more desirable for fast multi-person PE. However, the state-of-the-art
bottom-up methods are still computational intensive, usually requiring tens to
hundreds of GFLOPs. This is challenging for fast inference, especially on
resource-limited mobile devices.

DNN weight pruning He et al. (2019); Li et al. (2021b); Zhang et al. (2021),
as an effective model compression technique, has been widely adopted to remove
redundant weights in DNN models, reducing the required resources for both
model storage and computation, accelerating the model inference. The pruning
schemes mainly consist: 1) the flexible, unstructured pruning scheme that
prunes weights at arbitrary locations Han et al. (2016), 2) the regular,
structured pruning that prunes whole filters/channels for CONV layers He et
al. (2019), and 3) the compiler-assisted, fine-grained structured pruning
schemes Cai et al. (2020) that combining the high accuracy and high hardware
parallelism. Weight pruning provides a promising solution for accelerating
inference of multi-person PE on mobile devices while preserving good accuracy.

Model architecture is another critical factor that significantly affects PE
performance. For example, under similar computation FLOPs, the model depth has
a direct impact on model accuracy and inference latency. Moreover, the
unfriendly computation operators in a DNN model will introduce tremendous
execution overhead on resource-limited model devices but usually has a minor
impact on high-end server CPUs/GPUs. It is necessary to avoid using mobile-
unfriendly operators for on-mobile PE models.

Besides the DNN model inference, the data post-processing is another essential
part of the pose estimation process, which calculates and groups the keypoints
of each individual obtained from the DNN inference results. Different post-
processing algorithms will significantly affect the end-to-end running time
for the pose estimation. Thus, choosing appropriate post-processing algorithms
is also critical for both speed and accuracy consideration.

In this paper, to facilitate the fast and accurate multi-person PE on mobile
devices, we propose a framework incorporating architecture optimization and
weight pruning to find the desired model architecture and pruning
configurations to effectively accelerate inference of multi-person PE mobile
devices.

  * •

We propose an optimization framework for the multi-person PE acceleration. The
framework includes model architecture optimization and weight pruning to
facilitate fast and accurate end-to-end PE on mobile devices.

  * •

We propose several mobile-aware architecture optimization paradigms including
model depth determination, mobile-unfriendly operator replacement, and post-
processing selection.

  * •

We explore the sensitivity of different function blocks and different layer
types to the weight pruning and adopt suited pruning strategy for
optimization.

  * •

We achieve significant acceleration for multi-person PE inference on
smartphones while maintaining a high model accuracy.

  * •

Our achievements raise the promising future of various PE-related applications
on mobile devices.

##  2 Optimization Framework Design

The objective of our framework is to provide an optimization flow to
effectively accelerate PE models for fast inference on mobile devices while
maintaining competitive model accuracy. Our framework mainly contains two
parts, model architecture optimization, and weight pruning.

###  2.1 Model Architecture Optimization

####  2.1.1 Model Depth Determination

The models for multi-person PE tasks usually constitute several functional
blocks. We use the OpenPose Cao et al. (2019) model as an example. The
OpenPose model first uses a backbone model response for feature extraction.
The backbone depth will affect both final model accuracy and inference
latency. Generally, within a certain depth, a deeper backbone model provides a
better final accuracy, but leads to a higher inference latency. To avoid
significantly increasing the model computation FLOPs, we increase the number
of building blocks of the backbone model while shrinking the block width to
maintain similar computation FLOPs. As shown in Table 1, under similar
computation FLOPs, a deeper version of the ResNet50 model with 72 layers
provides a clear better final accuracy. But keep increasing model depth will
further improve the accuracy.

Besides, to determine the number of cascaded refinement blocks is also
important. The accuracy improvement obtained by adding refinement blocks
usually tends to be saturated after the first or second refinement block.
Thus, it is more economic to only add few refinement blocks, especially for an
inference speed-driven PE design on mobile devices.

####  2.1.2 Mobile-Unfriendly Operator Replacement

Using computation-friendly operators is often not the top concern when
designing networks pursuing high accuracy and using high-end server CPUs/GPUs.
But it is critical for the implementation on mobile devices. For example the
Swish activation function used in EfficientNet requires exponential
computation, which is mobile-unfriendly and significantly slows down the
inference speed. This even makes EfficientNet slower than ResNet50, which
requires higher computation FLOPs. Thus, our framework replace the unfriendly
operator with more friendly alternative (e.g., hard-tanh).

Some special kinds of convolution (CONV) layers are inefficient in mobile
inference such as the 7x7 CONV layer. In Lightweight OpenPose Osokin (2018),
they replace the 7x7 CONV layer with one 1x1 CONV layer followed by a 3x3 CONV
layer and a dilation 3x3 CONV layer. However, the dilation CONV layer requires
access to non-consecutive input data, which leads to bad locality and hence a
tremendous execution overhead. In our framework, we replace the 7x7 CONV layer
with three consecutive 3x3 CONV layers to preserve a similar receptive field
as the 7x7 CONV layer but with a much higher speed.

Model | Num of layers | GFLOPs | AP | mAP
---|---|---|---|---
ResNet$50$ | 24 | 13.68 | 0.41 | 0.682
ResNet$50_{2}$ | 36 | 16.05 | 0.423 | 0.69
ResNet$50_{3}$ | 72 | 16.09 | 0.436 | 0.693
ResNet$50_{4}$ | 147 | 16.46 | 0.433 | 0.69
Table 1: Accuracy under different depths of ResNet50 backbone.

####  2.1.3 Post-Processing Selection

After the DNN model inference, all PE methods require a post-processing step.
There are several different post-processing methods, and different methods
will result in significantly different post-processing times. Thus, selecting
an appropriate post-processing method is an essential problem for PE
implementations on mobile devices. To have a fair comparison of post-
processing time for different methods, we compared them using desktop CPU
since some methods are mobile-unfriendly and hard to be applied on mobile
devices.

In EfficientHRNet Christopher Neff (2020), it proposes to use the way same as
Associate Embedding, it needs to cluster the keypoints into multiple persons
according to tags got from the output of the network and match the heatmaps
for different resolutions, which is time-consuming. Although the
EfficientHRNet has smaller GFLOPs and higher accuracy, it takes around 2s for
post-processing on desktop CPU and dominates the overall end-to-end inference
time. As for HRNet Jingdong Wang and Xiao (2019), it used the keypoint
regression to improve keypoint detection and grouping, it still costs 0.2s. In
OpenPose Cao et al. (2019), it proposes to parse the heatmaps and PAFs
obtained from the network output by performing a set of bipartite matchings to
associate body part candidates and assemble them into full-body poses for all
people in the image. This post-processing takes about 0.02s. Both the first
two methods are not computational efficient and making the post-processing
become the bottleneck for the latency. The PAF post-processing used in
OpenPose is more desirable for the mobile-based implementation.

###  2.2 Weight Pruning

##### Pruning strategy:

In order to maintain accuracy while achieving a higher acceleration, we use
the fine-grained pruning schemes in our framework. Since different sizes and
types of DNN layers inherently have different redundancy, it is desirable to
adopt a layer-wise pruning ratio. Because sparse computation introduces
execution overhead, for one layer, only with a prune ratio that exceeds a
threshold can achieve actual acceleration. Our prune strategy for each layer
is either prune with a relatively high ratio or not prune. We established a
latency model to assist our prune ratio selection.

##### Sensitivity Analysis:

As different functional blocks (e.g., backbone, heatmap branch, PAF branch) in
a network exhibit different sensitivity to pruning, pruning on different
functional blocks lead to different impacts on the final accuracy. Therefore,
we conducted a sensitivity analysis on different functional blocks in the
model. We found that the backbone part has a higher tolerance for pruning and
can be pruned more aggressively. Compared with the heatmap branch, the PAF
branch is more sensitive to pruning, and the high pruning ratio will bring a
relatively higher accuracy drop, so we adopt a milder pruning strategy for the
PAF branch, that is, only prune a small number of layers or not prune.

##  3 Results and Demonstrations

Model Input Size Average Precision Num of Param. (M) FLOPs (G) Friendly Post-
Proce. PersonLab 1401 66.5 68.7 405.5 $\times$ Associate Embedding 512 65.5
138.9 222.3 $\times$ HRNet 512 64.4 28.5 38.9 $\times$ HigherHRNet 512 67.1
28.6 47.9 $\times$ EfficientHRNet $H_{-4}$ 384 35.7 3.7 2.1 $\times$
EfficientHRNet $H_{-3}$ 416 44.8 6.9 4.2 $\times$ OpenPose 368 48.6 52.3 136.1
$\surd$ Lightweight OpenPose 368 40 4.1 9 $\surd$ Ours ResNet$50_{3}$ (dense)
368 43.9 6.77 16.4 $\surd$ Ours ResNet$50_{3}$ (1.80$\times$) 368 42.0 3.76
9.7 $\surd$ Ours ResNet$50_{3}$ (2.15$\times$) 368 41.5 3.15 8.3 $\surd$ Ours
EfficientNet (2.01$\times$) 368 41.6 1.68 4.7 $\surd$

Table 2: Comparison with state-of-the-art pose estimation methods. Figure 1:
Demonstration of multi-person pose estimation on Samsung Galaxy S10
smartphone. Figure 2: Visual quality of multi-person pose estimation using our
optimized models. Figure 3: Model inference latency on Samsung Galaxy S10.

##### Experiment Setup:

We trained our models on MS COCO dataset Tsung-Yi Lin (2014) using PyTorch
API. The model inference latency is measured on a Samsung Galaxy S10
smartphone using a mobile GPU. Figure 1 shows the demonstrations of multi-
person pose estimation on the smartphone. We follow the design methodology of
OpenPose Cao et al. (2019) since the PAF the most mobile-friendly post-
processing compared to others. Our optimizations include model depth
determination, mobile-unfriendly operator replacement, and weight pruning.
Note that our optimization framework is general and can also be used to
optimize other design methodologies.

##### Accuracy & Visual Quality:

We compared our optimized models with representative PE methods Papandreou et
al. (2018); Alejandro Newell (2017); Jingdong Wang and Xiao (2019); Bowen
Cheng (2020); Christopher Neff (2020); Cao et al. (2019); Osokin (2018), as
shown in Table 2. Most of the methods either have huge computation FLOPs or
mobile-unfriendly post-processing, or both, making them not applicable to on-
mobile PE tasks. When compare with the Lightweight OpenPose (LWOP), which has
the potentiality to be deployed on mobile devices, our optimized models
achieve higher accuracy and lower FLOPs. Figure 2 shows the visual quality of
multi-person PE using our different optimized models. As shown in the figure,
we can maintain high visual quality with a significant reduction in
computation FLOPs.

##### Inference Latency:

Figure 3 shows the model inference latency comparison of our optimized models
and LWOP model on smartphones. As mentioned in Table 2 the LWOP model has
smaller computation FLOPs compared to our dense model and 1.80$\times$ pruned
model (i.e., 9G vs. 16.4G and 9.7G), but our optimized models are 1.42$\times$
and 2.07$\times$ faster, respectively. This is because we replaced the mobile-
unfriendly layers in the model. And our best model achieves 2.51$\times$
speedup with higher accuracy compared to LWOP.

##  4 Conclusion

We propose an optimization framework towards enabling fast and accurate multi-
person pose estimation on mobile devices. Our optimized models achieve both
higher accuracy and faster inference speed on mobile device.

## Acknowledgements

This project is partly supported by National Science Foundation (NSF) under
CNS-1739748, Army Research Office/Army Research Laboratory via grant
W911NF-20-1-0167 (YIP) to Northeastern University, a grant from Semiconductor
Research Corporation (SRC), and Jeffress Trust Awards in Interdisciplinary
Research.

## References

  * Alejandro Newell [2017] Jia Deng Alejandro Newell, Zhiao Huang.  Associative embedding: Endto-end learning for joint detection and grouping.  In Advances in Neural Information Processing Systems (NIPS), 2017\.
  * Bowen Cheng [2020] Jingdong Wang Honghui Shi Thomas S. Huang Lei Zhang Bowen Cheng, Bin Xiao.  Higherhrnet: Scale-aware representation learning for bottom-up human pose estimation.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * Cai et al. [2020] Yuxuan Cai, Hongjia Li, Geng Yuan, Wei Niu, Yanyu Li, Xulong Tang, Bin Ren, and Yanzhi Wang.  Yolobile: Real-time object detection on mobile devices via compression-compilation co-design.  arXiv preprint arXiv:2009.05697, 2020.
  * Cao et al. [2019] Z. Cao, G. Hidalgo Martinez, T. Simon, S. Wei, and Y. A. Sheikh.  Openpose: Realtime multi-person 2d pose estimation using part affinity fields.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2019\.
  * Christopher Neff [2020] Steven Furgurson Hamed Tabkhi Christopher Neff, Aneri Sheth.  Efficienthrnet: Efficient scaling for lightweight high-resolution multi-person pose estimation.  arXiv preprint arXiv:2007.08090, 2020.
  * Fang et al. [2017] Hao-Shu Fang, Shuqin Xie, Yu-Wing Tai, and Cewu Lu.  RMPE: Regional multi-person pose estimation.  In ICCV, 2017.
  * Han et al. [2016] Song Han, Huizi Mao, and William J. Dally.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  In International Conference on Learning Representations (ICLR), 2016\.
  * He et al. [2019] Yang He, Ping Liu, Ziwei Wang, Zhilan Hu, and Yi Yang.  Filter pruning via geometric median for deep convolutional neural networks acceleration.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 4340–4349, 2019.
  * Jingdong Wang and Xiao [2019] Tianheng Cheng Borui Jiang Chaorui Deng Yang Zhao Dong Liu Yadong Mu Mingkui Tan Xinggang Wang Wenyu Liu Jingdong Wang, Ke Sun and Bin Xiao.  Deep highresolution representation learning for human pose estimation.  IEEE transactions on pattern analysis and machine intelligence (T PAMI), 2019.
  * Li et al. [2021a] Hongjia Li, Geng Yuan, Wei Niu, Yuxuan Cai, Mengshu Sun, Zhengang Li, Bin Ren, Xue Lin, and Yanzhi Wang.  Real-time mobile acceleration of dnns: From computer vision to medical applications.  In 2021 26th Asia and South Pacific Design Automation Conference (ASP-DAC), pages 581–586. IEEE, 2021.
  * Li et al. [2021b] Zhengang Li, Geng Yuan, Wei Niu, Pu Zhao, Yanyu Li, Yuxuan Cai, Xuan Shen, Zheng Zhan, Zhenglun Kong, Qing Jin, Zhiyu Chen, Sijia Liu, Kaiyuan Yang, Bin Ren, Yanzhi Wang, and Xue Lin.  Npas: A compiler-aware framework of unified network pruning and architecture search for beyond real-time mobile acceleration.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2021.
  * Niu et al. [2020] Wei Niu, Zhenglun Kong, Geng Yuan, Weiwen Jiang, Jiexiong Guan, Caiwen Ding, Pu Zhao, Sijia Liu, Bin Ren, and Yanzhi Wang.  Achieving real-time execution of transformer-based large-scale models on mobile with compiler-aware neural architecture optimization.  arXiv preprint arXiv:2009.06823, 2020.
  * Osokin [2018] Daniil Osokin.  Real-time 2d multi-person pose estimation on cpu: Lightweight openpose.  In arXiv preprint arXiv:1811.12004, 2018.
  * Papandreou et al. [2018] George Papandreou, Tyler Zhu, Liang-Chieh Chen, Spyros Gidaris, Jonathan Tompson, and Kevin Murphy.  Personlab: Person pose estimation and instance segmentation with a bottom-up, part-based, geometric embedding model.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 269–286, 2018.
  * Tsung-Yi Lin [2014] Serge Belongie Lubomir Bourdev Ross Girshick James Hays Pietro Perona Deva Ramanan C. Lawrence Zitnick Piotr Dollár Tsung-Yi Lin, Michael Maire.  Microsoft coco: Common objects in context.  In Proceedings of the European Conference on Computer Vision (ECCV), 2014.
  * Zhang et al. [2021] Chengming Zhang, Geng Yuan, Wei Niu, Jiannan Tian, Sian Jin, Donglin Zhuang, Zhe Jiang, Yanzhi Wang, Bin Ren, Shuaiwen Leon Song, et al.  Clicktrain: Efficient and accurate end-to-end deep learning training via fine-grained architecture-preserving pruning.  In The 35th ACM International Conference on Supercomputing (ICS 2021), 2021.
  * Zhao et al. [2020] Pu Zhao, Wei Niu, Geng Yuan, Yuxuan Cai, Hsin-Hsuan Sung, Wujie Wen, Sijia Liu, Xipeng Shen, Bin Ren, Yanzhi Wang, et al.  Achieving real-time lidar 3d object detection on a mobile device.  arXiv preprint arXiv:2012.13801, 2020.