VARIATIONAL APPROACHES IN MOLECULAR ELECTROSTATICS, SURFACE
FORMATION, SHOCK CAPTURING AND NANO-TRANSISTORS
By
Langhua Hu

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
Applied Mathematics - Doctor of Philosophy
2013

ABSTRACT
VARIATIONAL APPROACHES IN MOLECULAR ELECTROSTATICS,
SURFACE FORMATION, SHOCK CAPTURING AND
NANO-TRANSISTORS
By
Langhua Hu
This dissertation covers several topics in Applied Mathematics, including nonlinear Poisson equation(NLPE) with application in electrostatics and solvation analysis for biological
system, partial diﬀerential equation(PDE) transform for hyperbolic conservation laws, high
order fractional PDE transform for molecular surface construction and Poisson-Kohn-Sham
equation for modeling geometric, thermal and tunneling eﬀects on nano-transistors.
Electrostatic interactions are ubiquitous in nature and fundamental for chemical, biological and material sciences. The Poisson equation is a widely accepted model for electrostatic
analysis. However, the Poisson equation is derived based on electric polarizations in a linear,
isotropic and homogeneous dielectric medium. We introduce a nonlinear Poisson equation
to take into consideration of hyperpolarization eﬀects due to intensive charges and possible
nonlinear, anisotropic and heterogeneous media. Variational principle is utilized to derive
the nonlinear Poisson model from an electrostatic energy functional. The proposed nonlinear
Poisson theory is extensively validated. A good agreement between our results and experimental data as well as theoretical results suggests that the proposed nonlinear Poisson model
is a potentially useful model for electrostatic analysis involving hyperpolarization eﬀects.
In our next work, we introduce the use of the PDE transform, paired with the Fourier
pseudospectral method (FPM), as a new approach for hyperbolic conservation law problems,
which remains an interesting and challenging task due to the diversity of physical origins

and complexity of the physical situations. A variety of standard benchmark test problems of
hyperbolic conservation laws is employed to systematically validate the performance of the
present PDE transform based FPM.
Furthermore, we study the high-order factional PDE transform based on fractional derivative with application in molecular surface generation. Fractional derivative or fractional calculus plays a signiﬁcant role in theoretical modeling of scientiﬁc and engineering problems.
However, only relatively low order fractional derivatives are used at present. Our work introduces arbitrarily high-order PDEs to describe fractional hyper-diﬀusions. The fractional
PDEs are constructed via fractional variational principle. The proposed high-order fractional
PDE transform are applied to the surface generation of proteins. Computational eﬃciency
of the present surface generation method is compared with the MSMS approach in Cartesian representation. Extensive numerical experiments and comparison with an established
surface model indicate that the proposed high-order fractional PDEs are robust, stable and
eﬃcient for biomolecular surface generation.
The last part of my work is in the ﬁled of nano-scale electronic transistors. The miniaturization of nano-scale electronic transistors, such as metal oxide semi- conductor ﬁeld eﬀect
transistors (MOSFETs), has given rise to a pressing demand in the new theoretical understanding and practical tactic for dealing with quantum mechanical eﬀects in integrated
circuits. We study the eﬀects of geometry of semiconductor-insulator interfaces, phononelectron interactions, and quantum tunneling of nano-transistors. Performances of nanotransistors are explored in terms of current-voltage (I-V) curves and quantized transport
energy proﬁles in a wide range of device parameters.

To my wife Chen, our kid Steven, parents, grandparents and extended family.

iv

ACKNOWLEDGMENTS

It is a long journey for me to ﬁnish my PhD dissertation. When I looked back over the process
of my PhD studies, I further realized how lucky I am since I have got a lot of valuable help
from many people. Within the limited space, I wish to acknowledge those who contribute to
my PhD journey. Before this, I would like to thank all my dissertation committee members,
Dr. Guowei Wei, Dr. Keith Promislow, Dr. Moxun Tang, Dr. Chichia Chiu and Dr. Yiying
Tong for taking their precious time to review and comment on my dissertation.
I’m indebted to my PhD advisor Dr. Guowei Wei, for accepting me as his student and
giving me consistent guidance and direction through my PhD studies. Dr. Wei has given
me step-by-step systematical training in programming and research. Not only he teaches
me knowledgeable in Mathematics, Biology, Physics and Chemistry but also trained me to
think sharply and critically in research. I appreciate him spending a lot of time designing
my projects, teaching me knowledge and skills and revising my manuscripts.
I’d like to thank former students in Dr.Wei’s group, Dr. Zhan Chen and Dr. Duan Chen,
for giving me many practical programming tools and research tips and kept helping me when
I met diﬃculties in researches. My thanks go to other Dr. Wei’s group alumni, Dr. Weihua
Geng, Dr. Sining Yu, Dr.Yuhui Sun, Dr. Yongcheng Zhou, and Dr. Shan Zhao for their
legacy code and research. Meanwhile, I thank Ms. Qiong Zheng, Dr. Siyang Yang, Dr.
Kelin Xia, Mr. Yin Cao, Ms. Weijuan Zhou and Ms. Jinkuyang Park for their help and
collaboration in my research and teaching.
I met generous faculty at MSU who gave me various kinds of training and help. To name
a few, I thank Dr. Baisheng Yan for the help in my ﬁrst year, Dr. Peiru Wu for the high
standard training in her class, Dr. Michael Feig for the training of Computational Biology
v

in his class and Dr. Andrew Christlieb for being my comprehensive exam member.
I am grateful to Ms. Barbara Miller, Graduate Secretary in the Department of Mathematics, for her generous assistance throughout my PhD study.
There are numerical colleagues and friends I should thank at MSU. To name a few, I
thank Mr. Shuzhuan Zheng, Mr. Xianfeng Ma, Mr. Li Zhan, Ms. Hong Qian, Mr. Xun
Wang, Mr. Lianzhang Bao, Mr. Junshan Lin, Mr. Jun Lai for their various help at diﬀerent
time.
Finally I would like to thank my family members, from the old generations of my grandparents and parents, to my generation and my son. Life was diﬃcult for the old generations
but they never gave up the hope and constantly made sacriﬁces for the next generation. I
especially owned my parents Aimin Hu and Qingying Zeng, my sister Fenghua Hu too much
for their long-time sacriﬁce and support. I utmost thanks go to my wife Chen Li, who gave
me constant support and encouragement and my one year old son Steven Li Hu, who gave
me a lot of joys and sense of responsibility to help me grow mature gradually as a father.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii

Chapter 1 Nonlinear Poisson equation for heterogeneous media
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Models and algorithms . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Electrostatic free energy functional . . . . . . . . . . . . .
1.2.2 Nonlinear Poisson equation . . . . . . . . . . . . . . . . .
1.2.3 Nonpolar energy . . . . . . . . . . . . . . . . . . . . . . .
1.2.4 Numerical Methods . . . . . . . . . . . . . . . . . . . . . .
1.2.4.1 Solution of nonlinear Poisson equation . . . . . .
1.2.4.2 Initialization of dielectric proﬁle ε . . . . . . . .
1.2.4.3 Evaluation of solvation free energy . . . . . . . .
1.2.4.4 Iteration procedure . . . . . . . . . . . . . . . . .
1.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Validation by one atom system . . . . . . . . . . . . . . .
1.3.2 Application to a set of 17 small molecules . . . . . . . . .
1.3.3 Application to a set of 20 proteins . . . . . . . . . . . . . .
1.3.4 Application to temperature eﬀects . . . . . . . . . . . . .
1.4 Chapter conclusion remarks . . . . . . . . . . . . . . . . . . . . .

1
1
4
4
8
11
15
15
17
17
18
19
20
23
28
31
33

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

Chapter 2 PDE transform for hyperbolic conservation laws . . . . . . . . .
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Theory and algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Arbitrarily high order nonlinear PDEs . . . . . . . . . . . . . . . . .
2.2.2 PDE transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Numerical algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Numerical test and validation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Scalar conservation law systems . . . . . . . . . . . . . . . . . . . . .
2.3.1.1 Example 1 (Linear advection equation with Sine-Gaussian
wavepacket) . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1.2 Example 2 (Linear advection equation with wave combination)
2.3.1.3 Example 3 (Inviscid Burgers’ equation) . . . . . . . . . . .
2.3.1.4 Example 4 (Non-convex ﬂux) . . . . . . . . . . . . . . . . .
2.3.2 1D Euler systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2.1 Example 5 (Sod’s and Lax’s problems) . . . . . . . . . . .
2.3.2.2 Example 6 (1D Shock-entropy interaction) . . . . . . . . .
2.3.2.3 Example 7 (Shu-Osher’s problem) . . . . . . . . . . . . . .

vii

35
35
42
42
46
48
51
53
54
58
62
64
67
68
71
77

2.3.3

2.4

2D Euler systems . . . .
2.3.3.1 Example 8 (2D
2.3.3.2 Example 9 (2D
Chapter conclusion remarks . .

. . . . . . . . . . . . . . . .
Shock-entropy interaction)
Shock-vortex interaction) .
. . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Chapter 3
3.1
3.2

3.3

3.4

3.5

High order fractional PDE transform for molecular surface
construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Theory and algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Fractional derivatives . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Fractional ﬁnite diﬀerence schemes . . . . . . . . . . . . . . . . . . .
3.2.2.1 Fractional Fourier schemes . . . . . . . . . . . . . . . . . .
3.2.2.2 Integral forms . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 High-order fractional PDEs and fractional PDE transform . . . . . .
3.2.3.1 Variational derivation of high-order fractional PDEs . . . .
3.2.3.2 Fractional hyperdiﬀusions derivation of high-order fractional
PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 Numerical algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . .
Numerical test and validation . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Initial data for evolution PDEs . . . . . . . . . . . . . . . . . . . . .
3.3.1.1 Type I: Piecewise-constant initial values . . . . . . . . . . .
3.3.1.2 Type II: Maximum Gaussian initial values . . . . . . . . . .
3.3.1.3 Type III: Summation Gaussian initial values . . . . . . . . .
3.3.2 Eﬀects of fractional order and propagation time . . . . . . . . . . . .
3.3.2.1 Test on a three-atom system . . . . . . . . . . . . . . . . .
3.3.2.2 Test on a four-atom system . . . . . . . . . . . . . . . . . .
3.3.2.3 Test on two proteins . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Computational eﬃciency . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Surface areas and surface enclosed volumes . . . . . . . . . . . . . . .
Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Electrostatic analysis of 15 proteins . . . . . . . . . . . . . . . . . .
3.4.1.1 Electrostatic potential . . . . . . . . . . . . . . . . . . . . .
3.4.1.2 Electrostatic solvation free energy . . . . . . . . . . . . . .
3.4.2 Solvation analysis of 17 small compounds . . . . . . . . . . . . . . .
Chapter conclusion remarks . . . . . . . . . . . . . . . . . . . . . . . . . . .

79
80
83
87

100
101
104
104
104
105
105
105
108
111
114
115
116
119
119
119
121
121
124

Chapter 4
4.1
4.2

Modeling of Impact of geometric, thermal and tunneling eﬀects
on nano-transistors . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Theoretical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Free energy of electrons . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1.1 Kinetic energy . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1.2 Potential energy . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Electrostatic free energy of the system . . . . . . . . . . . . . . . . .
4.2.3 Total free energy functional . . . . . . . . . . . . . . . . . . . . . . .

91
91
92
92
93
94
96
98
98

viii

127
127
131
131
131
132
134
135

4.3

4.4

4.5

4.2.4 The Poisson equation . . . . . . . . . . . . . . . . . . .
4.2.5 Generalized Kohn-Sham equation . . . . . . . . . . . .
Computational algorithms . . . . . . . . . . . . . . . . . . . .
4.3.1 The Poisson equation . . . . . . . . . . . . . . . . . . .
4.3.2 Scattering problem . . . . . . . . . . . . . . . . . . . .
4.3.3 Modeling the phonon-electron interaction . . . . . . . .
4.3.4 Self-consistent iteration . . . . . . . . . . . . . . . . . .
Numerical experiments and discussions . . . . . . . . . . . . .
4.4.1 A four-gate MOSFET model . . . . . . . . . . . . . . .
4.4.2 General results of electric ﬁeld and electron density . .
4.4.3 Performances of nano-MOSFETs with diﬀerent channel
4.4.4 Quantum tunneling eﬀects . . . . . . . . . . . . . . . .
4.4.5 Phonon-electron interactions . . . . . . . . . . . . . . .
Chapter conclusion remarks . . . . . . . . . . . . . . . . . . .

. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
geometries
. . . . . . .
. . . . . . .
. . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

135
136
137
137
139
141
142
144
144
148
151
153
158
160

Chapter 5 Thesis achievements and future work . . . . . . . . . . . . . . . . 163
5.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
5.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
BIBLIOGRAPHY

. . . . . . . . . . . . .

ix

. . . . . . . . . . . . . . . . . 170

LIST OF TABLES

Table 1.1

The convergence of ∆Gp with respect to grid size h for one-atom
system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

Table 1.2

Comparison of atomic radii and adjusted values (˚). . . . . . . . . .
A

23

Table 1.3

Comparison of free energies (kcal/mol) for 17 compounds between
nonlinear model and experimental data. . . . . . . . . . . . . . . . .

27

Errors of solvation free energies for a set of 17 small compounds obtained from 3 nonlinear models based on adjusted radius in Table
1.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

Electrostatic solvation free energies (kcal/mol) for 20 proteins by diﬀerent models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Comparison of RMSE and ∆∆G predicted by the classical Poisson
model and the NLPE model for 21 compounds (The average ∆∆G of
experimental data is 1.13 kcal/mol). . . . . . . . . . . . . . . . . . .

33

Temperature eﬀect for 21 compounds between the NLPE model and
experimental results. . . . . . . . . . . . . . . . . . . . . . . . . . .

34

The order and propagation time of the PDE transform used in test
examples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

Table 2.2

L∞ error for long-time integration of the Sine-Gaussian wavepacket.

57

Table 3.1

Comparison of CPU costs (in seconds) for surface generation in Cartesian meshes by the fractional PDE transform (FPDE) with α = 11.5,
t = 104 and the MSMS approach. Na : total number of atoms; Ts :
CPU time for MSMS to generate surface in 2D Lagrangian representation; Tc : CPU time for converting the Lagrangian representation
to 3D Eulerian representation. Ttot : Total CPU cost for FPDE and
MSMS; h: grid resolution. . . . . . . . . . . . . . . . . . . . . . . . 117

Table 1.4

Table 1.5

Table 1.6

Table 1.7

Table 2.1

x

Table 3.2

˚
Comparison of surface areas (unit: A2 ), surface enclosed volumes
3 ) and electrostatic free energies of solvation ∆G (unit:
(unit: ˚
A
p
kcal/mol) with surfaces generated by the fractional PDE transform
(FPDE) and the MSMS approach. The grid resolution is 0.5 ˚ . . . 118
A

Table 3.3

Computation of solvation free energies ∆G (kcal/mol) for the 17 small compounds based on the molecular surface generated by fractional
PDE transform with order α = 11.5 and propagation time t = 104 .
Gnp : nonpolar component; ∆Gp : polar component; Exptl: experimental data; Error: ∆Gp - Exptl. . . . . . . . . . . . . . . . . . . . . 122

Table 4.1

Tunneling ratio for the cylinder channel at diﬀerent VGate (VDS =0.4V).157

xi

LIST OF FIGURES

Figure 1.1

The cross section proﬁles of ε(| φ|) along the y-axis at x = 0 and
z = 0 at diﬀerent α values for one-atom system. For interpretation
of the references to color in this and all other ﬁgures, the reader is
referred to the electronic version of this dissertation. . . . . . . . .

10

The cross line proﬁles of S and 1 − S along the y-axis at x = 0 and
z = 0 for one-atom system. . . . . . . . . . . . . . . . . . . . . . .

12

Figure 1.3

The ﬂowchart for the solution of the nonlinear Poisson equation. . .

19

Figure 1.4

The convergence of electrostatic solvation free energies ∆Gp over
number of iteration (N) for one-atom system. . . . . . . . . . . . .

20

The volume of the nonlinear Poisson model at a wide range of α values
for one-atom system with atomic radius 1 ˚, εm = 1 and εs = 80. .
A

21

The electrostatic solvation energy of the nonlinear Poisson model at
a wide range of α values for one-atom system with atomic radius 1
˚, εm = 1 and εs = 80. . . . . . . . . . . . . . . . . . . . . . . . .
A

22

The rescaled atomic radii for the one-atom system with ﬁxed ε = εm
within the solute domain. . . . . . . . . . . . . . . . . . . . . . . .

24

Correlation of total solvation free energies of 17 compounds between
the nonlinear dielectric model and diﬀerent data. (a) Nonlinear dielectric model and experimental data; (b) Nonlinear dielectric model
and DGSM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

Correlation between the results obtained by the MIBPB-III [67] and
the nonlinear dielectric model for electrostatic solvation free energies
of 20 proteins. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

Surface electrostatic potentials of four molecules at their isosurface
value S = 0.5. (a) small compound: diethyl propanedioate; (b)small
compound: n,n-dimethyl-p-methoxybenzamide; (c)protein: 1vjw; (d)
protein: 2pde. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Figure 1.2

Figure 1.5

Figure 1.6

Figure 1.7

Figure 1.8

Figure 1.9

Figure 1.10

xii

Figure 1.11

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

Figure 2.8

Figure 2.9

Figure 2.10

Correlation of solvation free energies predicted by the NLPE model
and experimental data for 21 compounds. . . . . . . . . . . . . . .

30

Results for κ = 5, 10 from the PDE transform for the advection equation with the Sine-Gaussian wave packet (t=100, ∆t = 10−4 ). .

55

Results for κ = 20 from the PDE transform for the advection equation
with the Sine-Gaussian wave packet (t=100, ∆t = 10−4 ). . . . . . .

56

Results from the PDE transform of various orders for the advection
equation ( t = 8, ∆t = 0.001, 256 grid points). . . . . . . . . . . . .

60

Results from the 6th-order PDE transform for the inviscid Burgers’
equation at time t = 2. (a) ul = 1, ur = 0, 257 grid points; (b)
ul = 0, ur = 1, 129 grid points; . . . . . . . . . . . . . . . . . . . .

61

Comparison between numerical results from the PDE transform and
FPM-RSK method for the problem with a non-convex ﬂux ( 129 grid
points, t = 0.04, ∆t = 0.0005). . . . . . . . . . . . . . . . . . . . . .

65

Comparison between numerical results from the PDE transform and
FPM-RSK method for the problem with a non-convex ﬂux ( 65 grid
points, t = 0.04, ∆t = 0.0005). . . . . . . . . . . . . . . . . . . . . .

66

Numerical results from the 6th-order PDE transform for Sod’s problem (t = 1.5, ∆t = 0.02, 129 grid points). (a) Density; (b) Pressure.

69

Comparison of numerical results from the 6th-order PDE transform
and the FPM-RSK for Lax’s problem (t = 1.5, ∆t = 0.02, 129 grid
points). (a) Density from the PDE transform; (b) Pressure from the
PDE transform; (c) Density from the FPM-RSK; (d) Pressure from
the FPM-RSK. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

The amplitude of entropy waves in the 1D shock-entropy interaction
problem for diﬀerent frequencies resulting from the 12th-order PDE
transform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

The amplitude of entropy waves in the 1D shock-entropy interaction problem for κ = 32 resulting from PDE transforms of low and
extremely high-order. (a) The 6th-order PDE transform; (b) The
40th-order PDE transform. . . . . . . . . . . . . . . . . . . . . . . .

73

xiii

Figure 2.11

Figure 2.12

Figure 2.13

Figure 2.14

Solution to the Shu-Osher problem (t = 0.47, ∆t = 0.001, 129 grid
points). (a) Obtained with the 2nd-order PDE transform; (b) Obtained with the 10th-order PDE transform. . . . . . . . . . . . . .

77

The amplitude of entropy waves in the 2D shock-entropy interaction
problem obtained from PDE transforms of diﬀerent orders. . . . . .

81

The pressure proﬁle of 2D shock-vortex interaction problem from the
10th-order PDE transform at t = 0.05, 0.2 (20 contours). . . . . . . .

84

The pressure proﬁle of 2D shock-vortex interaction problem from the
10th-order PDE transform at t = 0.35, 0.6 (20 contours). . . . . . . .

85

Figure 3.1

Frequency response of Eq. (3.34) with t = 102 and diﬀerent orders.
(a) Over all amplitude; (b) Real component; (c) Imaginary component. Red: α = 1.5; blue: α = 5, 5; and purple: α = 11.5. . . . . . . 107

Figure 3.2

Real component of frequency response of Eq. (3.34) in 1D, with
diﬀerent values of α and t: (a) α = 11.2; (b) α = 11.5; (b) α = 12.0.
Red: t = 106 ; blue: t = 104 ; and purple: t = 102 . . . . . . . . . . . 109

Figure 3.3

Real component of Fourier response of Eq. (3.34) in 2D for α = 11.5
and t = 102 . (a) function value view; (b) contour view. . . . . . . . 110

Figure 3.4

Isosurfaces of a three-atom system after the fractional PDE transform
of diﬀerent fractional orders with the piecewise constant initial value.
The isovalue is taken as u(r, t) = 1 and propagation time is t = 102 .
(a) Initial value; (b) α = 1.5; (c) α = 5.5; (d) α = 11.5. . . . . . . . 111

Figure 3.5

Isosurface of a four-atom system after fractional PDE transform of
order α = 11.5 and diﬀerent propagation time, with the piecewise
constant initial value. The isovalue is taken as u(r, t) = 1 (a) initial
value; (b) t = 1; (c) t = 104 . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 3.6

Isosurfaces of a four-atom system after fractional PDE transform of
order α = 11.5 and diﬀerent propagation time, with maximum Gaussian initial value. The isovalue is taken as u(r, t) = 1. (a) Initial
value; (b) t = 1; (c) t = 104 . . . . . . . . . . . . . . . . . . . . . . . 113

xiv

Figure 3.7

Isosurfaces of a four-atom system after fractional PDE transform of
oder α = 11.5 and diﬀerent propagation time, with the summation
Gaussian initial value. The isovalue is taken as u(r, t) = 1 for (a)-(c)
and u(r, t) = 2 for (d)-(f). (a) Initiale value; (b) t = 1; (c) t = 104 ;
(d) initial value; (e) t = 1; (f) t = 104 . . . . . . . . . . . . . . . . . 113

Figure 3.8

Isosurface of a four-atom system with initial value of Type III (Summation Gaussian initial values) at diﬀerent κ. The isovalue is taken
as u(r, 0) = 1. (a) κ = 1; (b) κ = 1/9. . . . . . . . . . . . . . . . . 114

Figure 3.9

Surface of a protein (PDB code: 1R69) after fractional PDE transform of propagation time t = 104 and diﬀerent orders. The isovalue
is taken as u(r, t) = 1. (a) Initial value (Type I); (b) α = 1.5; (c)
α = 11.5; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Figure 3.10

Surface of a protein (PDB code: 1FCA) after fractional PDE transform of order α = 11.5 at diﬀerent propagation time. The isovalue is
taken as u(r, t) = 1. (a) initial value (Type I); (b) t = 100 ; (c) t = 104 ; 115

Figure 3.11

Surface electrostatic potentials of three proteins based on the molecular surface after fractional PDE transform with α = 11.5 and t = 104 .
(a) PDB code: 1VII; (b) PDB code: 1FCA; (3) PDB code:1R69. Red:
negative potential; blue: positive potential. . . . . . . . . . . . . . . 120

Figure 3.12

Illustration of surface electrostatic potentials of nine small compounds at their corresponding isosurfaces u(r, 104 ) = 0.5. (a) benzyl bromide; (b) bis-2-chloroethyl ether; (c) 1,1-diacetoxyethane; (d) diethyl
propanedioate; (e) diethyl sulﬁde; (f) dimethoxymethane; (g) m-bis
benzene; (h) N,N-4-trimethylbenzamide; (i) phenyl formate. . . . . . 126

Figure 4.1

The ﬂowchart of Gummel-like self-consistent iteration. . . . . . . . . 144

Figure 4.2

Illustration of a four-gate MOSFET with a cylinder channel: (a)
Overall 3D conﬁguration; (b) View of the x-z cross section. . . . . . 145

Figure 4.3

(a) y-z cross-section of a channel: circle; (b) y-z cross-section of channel: ﬂower-like shape. . . . . . . . . . . . . . . . . . . . . . . . . . 146

Figure 4.4

Demonstration of the simpliﬁed computational domain for the MOSFET. Dark: insulator layer; grey: doped channel region; and white:
un-doped channel region. Gate voltages are applied at boundary
segments EF and IJ, and source-drain voltages are applied on the
boundary segments GH and LK. . . . . . . . . . . . . . . . . . . . . 147
xv

Figure 4.5

Electrostatic and electron density proﬁles from the x-z cross-section
view at y = 2.5nm in the silicon area of the MOSFET with the
cylinder channel (Source-drain voltage VDS = 0.4V and gate voltage
VGate = 0.4V). (a) Electrostatic potential energy (eV); (b) Electron
density (1e26/m3 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

Figure 4.6

Subband energy proﬁle along the transport direction under gate voltages VGate = 0.1V(Dashed line) and VGate = 0.4V(Solid line) for the
MOSFET with the cylinder channel (source-drain voltages is ﬁxed at
VDS = 0.4V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

Figure 4.7

Current-voltage (I - V) curves of the MOSFET with the cylinder
channel. (a) I-VDS curve (gate voltage ﬁxed at VGate =0.4V); (b)
I-VGate curve (source-drain voltage ﬁxed at VDS =0.4V). . . . . . . 150

Figure 4.8

I-VGate curves of MOSFETs with diﬀerent geometries of channel
(VDS =0.4V). (a) Four types of cross-sections channels; (b) Corresponding I-VGate curves. . . . . . . . . . . . . . . . . . . . . . . . . 151

Figure 4.9

The I-VGate curves of MOSFETs with diﬀerent channel cross-sections
but with the same cross-section area ( source-drain voltage ﬁxed at
VDS =0.4V). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

Figure 4.10

Subband energy proﬁles along the transport direction for MOSFETs
with diﬀerent geometries of channel cross-sections(source-drain voltages is ﬁxed at VDS = 0.4V and gate voltages VGate = 0.4V). (a)
First subband energies of diﬀerent channels; (b) Diﬀerence of ﬁrst
two subband energies of each channel. . . . . . . . . . . . . . . . . . 154

Figure 4.11

Relation between electron current (µA) and electron energy E (eV) in
the MOSFET with the cylinder channel. (a)Gate voltage and sourcedrain voltage are respectively VGate = 0.0V and VDS = 0.4V. The red
dot indicates the current w.hen electron energies equal to the maximum energy barrier Emax (horizontal coordinate). (b)under various
gate voltages VGate (source-drain voltages ﬁxed at VDS = 0.4V). The
red dots from right to left correspond to VGate = 0, 0.05, 0.1, ..., 0.4V. 156

Figure 4.12

Relation of electrons tunneling ratio and gate voltages VGate for nanoMOSFETs with cylinder and ﬂower-like channels (source-drain voltage ﬁxed at VDS = 0.4V). . . . . . . . . . . . . . . . . . . . . . . . . 158

xvi

Figure 4.13

Eﬀects of phonon-electron interactions in I-V curves and electron tunneling ratio for the cylinder channel. (a) I-VGate curves with diﬀerent electron-phonon interactions; (b) Relations of electron tunneling
ratios and gate voltages VGate with diﬀerent electron-phonon interactions (source-drain volate ﬁxed at VDS =0.4V). . . . . . . . . . . . 159

Figure 4.14

Temperature eﬀects on electron currents under diﬀerent phonon-electron
interactions for the cylinder channel (VDS =0.4V, VGate =0.3V). (a)
Current-temperature curves; (b) Electron tunneling ratios under different temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . 160

xvii

Chapter 1
Nonlinear Poisson equation for
heterogeneous media

1.1

Introduction

Electrostatic interactions are ubiquitous in nature and fundamental for chemical, biological
and material sciences. Among the various components of molecular interactions, electrostatic interactions are of special importance [174, 46, 53, 83, 146, 173, 147, 61, 148, 10, 58, 11]
because of their long range and inﬂuence on polar or charged molecules – including water,
aqueous ions, proteins, nucleic acids, lipid bilayers, sugars, etc. For example, proteins are
made up of 20 types of amino acids, 11 of which are charged or polar in neutral solution.
Nucleic acids contain long stretches of negative charges from the phosphate groups in nucleotides. In addition, the double-strained DNA chains are found to have an electron charge
density as high as one charge per 1.7 ˚. Finally, sugars and related glycosaminoglycans can
A
possess some of the highest charge densities of any biomolecules because of the presence
of numerous negative functionalities, including carboxylate and sulfate groups. Therefore,
electrostatic interactions are some of the most important aspects that determine the physical and chemical properties of biomolecules, such as protein folding, protein-DNA binding,
transcription, translation, gene expression, and regulation, etc. As most biological processes occur in aquatic environments, electrostatic solute-solvent interactions are of paramount
1

importance in the exploration of biological mechanism, the analysis of macromolecular behavior, and the modeling of the intramolecular and inter- molecular interactions of biological
complexes.
A wide variety of theoretical approaches, ranging from quantum mechanical ab initio
methods, classic Maxwell’s theory of electromagnetism, generalized Born algorithms, to
phenomenological modiﬁcations of Coulomb’s law, has been used for electrostatic analysis. Quantum mechanical methods are generally accurate, but are too expensive to be used
for large chemical and biological systems. Generalized Born algorithms are fast, but depend
on other methods for calibrations. The Poisson equation, derived by using the Gauss’s law
and linear polarization, provides a relatively simple and accurate, while much less expensive
description of electrostatic interactions for a given charge source density. It works quite well
for moderately charged small molecules such as single amino acids [130].
One of most important applications of the Poisson equation is the solvation analysis
in the framework of implicit solvent models [175, 92, 139, 37]. Solvation involves interactions between solute molecules and solvent molecules or ions in the aqueous environment.
Since 65-90% of cellular mass is water, all important biological processes, such as protein
ligand binding, ion transport, signal transduction, gene regulation, transcription, and translation, occur in aqueous environments under physiological conditions. Therefore, one cannot
overemphasize the importance of solvation analysis to chemical and biological sciences. The
Poisson equation, or the Poisson-Boltzmann equation when the salt is present, is a main
workhorse in the electrostatic free energy estimation of solvation because of its balance of
model accuracy and computational cost.
However, the Poisson equation is strictly valid for the linear polarization in an isotropic
and homogeneous dielectric medium. In contrast, one frequently encounters heavily charged,
2

nonlinear, anisotropic and inhomogeneous materials and strong electric ﬁelds in real world
problems, for which hyperpolarizations become signiﬁcant. For example, hyperpolarizability
often manifests itself as the Pockels eﬀect and/or the Kerr eﬀect [25, 87, 72, 117, 122].
Nonlinear dielectric response is important for highly charged molecules such as DNA or in
protein active sites. It has also been found in composite ferroelectric materials [137, 108, 156].
Nevertheless, for most complex material and devices, such as nanoﬂuidics, fuel and solar cells,
electronic transistors, ion channels, molecular motors, etc., the hyperpolarizability is usually
neglected in their electrostatic analysis. Therefore, there is a need to consider the nonlinear
modiﬁcation of the Poisson model for practical applications.
In the solvation analysis using the Poisson model or the Poisson-Boltzmann model in
the presence of salt, the solute is described in molecular or atomic detail, while the solvent
is described by a dielectric continuum, or a mean ﬁeld approximation. The media in this
problem are described by variable coeﬃcients and regarded as heterogeneous. The solventsolute boundary is often prescribed by solvent accessible, solvent excluded or van der Waals
surfaces, which gives rise to an unphysical sharp jump in the dielectric proﬁle. Various
modiﬁcations of the Poisson-Boltzmann model have been considered in the literature, including dipolar solvents [7, 1, 63], multipole eﬀects [145, 102], multiple dielectric constants
and multivalent ions [138], solvent-solute interaction [106, 37, 55], and steric eﬀects [24].
Warshel presented a protein-dipole Langevin-dipole (PDLD) model [172], and Freed oﬀered
a Langevin-Debye theory [73].
The emphasis in most earlier work is on the solvent, while nonlinear eﬀects in the vicinity
of the solute boundary, and their impact to heterogeneous and/or anisotropic media have
received relatively less attention. The present work introduces a nonlinear Poisson equation
to model nonlinear eﬀects in solvent-solute systems in the absence of salt.
3

1.2
1.2.1

Models and algorithms
Electrostatic free energy functional

Biomolecules are typically dielectric materials which do not conduct electric current, while
their electrons redistributes in response to an applied electric ﬁeld. The charge redistribution
leads to so called dielectric polarization. The polarization density P is proportional to the
electric ﬁeld E and can be measured by susceptibility χe

P = ε0 χe E,

(1.1)

where ε0 is the electric permittivity of free space. The susceptibility of a medium is related
to its relative permittivity ε by χe = ε−1, which is set to 0 at vacuum where there is nothing
to be polarized. The electric displacement D is related to the polarization density P by

D = ε0 E + P = ε0 (1 + χe )E = ε0 εE.

(1.2)

By Gauss’s law, one has the Poisson equation

·D=

· ε0 εE = −

· ε0 ε φ(r) = ρ(r),

(1.3)

where φ(r) is the electrostatic potential and ρ(r) is the charge density. The relation of
E = − φ(r) has been used.
An essential assumption in the above derivation of the Poisson equation is that the
polarization density is linearly proportional to the electric ﬁeld E. However, this is not true
for highly charged nonlinear materials under a strong electric ﬁeld or in an inhomogeneous
4

medium with a complex interface geometry [87, 72]. Note that according to Coulomb’s law,
the electric ﬁeld diverges in the close vicinity of a point charge. In general, the polarization
density is of the form

P = ε0 χ(1)

(1)

E + ε0 χ(2)

+ ε0 χ(3)

(3)

(2)

EEE + · · · ,

EE
(1.4)

where χ(n) is an nth order susceptibility and is an (n + 1)th order hyperpolarizability tensor
and

(n)

is an nth order inner product for appropriate tensorial contraction. The second

order susceptibility describes the second order harmonic generation and Pockels eﬀect, while
the third order susceptibility represents the Kerr eﬀect [25]. As a vector in R3 , the components of P are given by

(1)

(2)

Pα = ε0 χαβ Eβ + ε0 χαβγ Eβ Eγ
(3)

+ ε0 χαβγδ Eβ Eγ Eδ + · · · ,

(1.5)

where the Einstein summation convention is used, i.e., repeated indices are summed over
all of its possible values. A major diﬃculty in the application of Eq. (1.5) is that there
are too many tensorial components in dealing with high order tensors. Symmetry of the
material can be utilized to reduce the number of independent parameters. However, for
complex molecules, such as proteins, their symmetry is normally very low and there is little
simpliﬁcation.
The polarization density described in Eq. (1.5) suggests that the linear Poisson equation
is a ﬁrst order approximation. In general, the Poisson equation should be nonlinear in E or

5

φ(r) so as to partially account for the eﬀect of hyperpolarizations in solvation analysis.
To derive a nonlinear Poisson equation, we consider an energy functional of the form

Gp =

(Λ ( φ(r)) + g(φ(r)) dr,

(1.6)

where g(φ(r)) is a source term and Λ(·) is an appropriate function to represent nonlinear
material properties. Related free energy functionals for electrostatic analysis were discussed
in the literature [146, 71, 180]. Energy functionals and variational approaches are now being
widely used to derive new and more sophisticated models of solvation [180, 37, 38, 39].
In general, one may set Λ(·) as a polynomial of the form

Λ( φ(r)) = aj,k | k φ(r)|j ,

j, k = 1, 2, · · · .

(1.7)

In this expression, the high order electric ﬁeld eﬀects are included. However, similar to Eq.
(1.5), to use Eq. (1.7), one faces the diﬃculty of determining parameters aj,k . The lack of
detailed information about chemical and physical properties of a molecule in solvent limits
the usage of Eq. (1.7).
In the present work, we consider part of the nonlinear polarization eﬀect in the solvation
analysis. To avoid the involvement of susceptibility parameters and maintain the simplicity
of the implicit solvent model, and to still account for the major hyperpolarizability impact
to the solvation, we propose a simple form

εm
| φ(r)|2
2
kB T
α
+
(εs − εm ) ln 1 +
| φ(r)|2 ,
α
2kB T

Λ(| φ(r)|) =

6

(1.8)

and

g(φ(r)) = −ρm (r)φ(r)

(1.9)

where εs and εm are the permittivities, i.e. dielectric constants, of the solvent and solute
respectively, kB is the Boltzmann constant, T is the temperature, α is a parameter to adjust
the strength of the hyperpolarizability, and ρm (r) is the molecular charge density
Nm

Qi δ(r − ri ),

ρm (r) =

(1.10)

i

with Qi being the partial charge on an atom located at position ri and Nm is the total
number of atoms in the molecule. The ﬁrst term in Eq. (2.11) and the choice of g(φ(r))
in Eq. (1.9) are standard for the classic electrostatic theory. Whereas, the second term in
Eq. (2.11) represents the partial eﬀect of the hyperpolarizability near the solvent-solute
boundary.
Certainly, the free energy functional given in Eq. (2.11) is not unique. For example, two
alternative forms are given by

εm
| φ(r)|2
2
1
k T
− B (εs − εm )
α | φ(r)|2 − 1 ,
α
1 + 2k T
B

Λ(| φ(r)|) =

(1.11)

and
| φ(r)|2
−α
2kB T
(εs − εm ) e



Λ(| φ(r)|) =

εm
k T
| φ(r)|2 − B
2
α

7


− 1 .

(1.12)

These energy functionals describe the α and temperature dependence of the hyperpolarizability. The impact of the hyperpolarization decays as the temperature increases [117, 122].

1.2.2

Nonlinear Poisson equation

Due to the simple form in Eqs. (2.11), (1.11) and (1.12), the energy functional in Eq. (2.11)
can be optimized with respect to the electrostatic potential φ by using the Euler-Lagrange
equation,

−

·

∂g
∂Λ
+
= 0.
∂ φ ∂φ

(1.13)

Based on the energy functionals in Eqs. (2.11) and (1.11) and Euler-Lagrange equation
(1.13), we obtain the nonlinear Poisson equation (NLPE)

−

· (ε(| φ(r)|) φ(r)) = ρm (r),

(1.14)

where the nonlinear dielectric function is given by

ε(| φ(r)|) = εm +

εs − εm
| φ(r)|2
1 + α 2k T
B

n

(1.15)

with n = 1 and 2 for energy functionals in Eqs. (2.11) and (1.11), respectively. For Eq.
(1.12), the nonlinear dielectric function is given by

−α

ε(| φ(r)|) = εm + (εs − εm ) e

8

| φ(r)|2
2kB T .

(1.16)

It has been shown in the literature [122] that the hyperpolarization decays exponentially as
the temperature increases.
The NLPE (1.14) has the same structure as the classic Poisson equation (1.2). However,
the dielectric functions in Eqs. (1.15) and (1.16) are no longer discontinuous. To understand
the behavior of ε(| φ|), we ﬁrst analyze two extreme cases,

lim ε(| φ|)→εs

(1.17)

lim

(1.18)

| φ|→0

| φ|→∞

ε(| φ|)→εm .

Obviously, Eqs. (1.17) and (1.18) give rise to desirable asymptotic dielectric values in the
solvent and the solute domains, respectively. It is easy to see that εm < ε(| φ|) < εs for
intermediate | φ| values. In fact, the dielectric functions ε(φ) in Eqs. (1.17) and (1.18)
undergo a continuous transition from εm to εs at the solvent-solute interface when | φ| is
continuous. An interesting issue is the transition length scale of ε(| φ|). By Coulomb’s law,
the electrostatic potential φ resulting from a point charge, supposed at the origin, decays as
1/r. Therefore, the inﬂuence of | φ|2 decays as 1/r4 , which is much faster than the decay
of the electrostatic potential and leads to a relatively short transition region of the dielectric
function.
Figure 1.1 illustrates the smooth transition of ε(| φ|) from εs to εm . Here, εs and εm
are set to their conventional values, i.e., 80 and 1, respectively. To understand the impact of
the hyperpolarizability, we plot ε(| φ|) with respect to a few α values in Figure 1.1. Larger
α values lead to a slower transition of ε(| φ|) to its asymptotic value εs . When there is
little hyperpolarization, i.e., α → 0, ε(| φ|) undergoes a sharp jump from its molecular
value εm to its solvent value εs . We therefore recover the classic Poisson equation model for
9

120
α=0.1
α=1
α=10
α=50

100

80

ε

60

40

20

0
−6

−4

−2
0
2
Y (angstrom)

4

6

Figure 1.1: The cross section proﬁles of ε(| φ|) along the y-axis at x = 0 and z = 0 at
diﬀerent α values for one-atom system. For interpretation of the references to color in this
and all other ﬁgures, the reader is referred to the electronic version of this dissertation.
the solvation analysis. However, the classic Poisson equation is ill-posed unless appropriate
interface jump conditions are enforced [67]. Mathematically, the classic Poisson equation is
very diﬀerent from the proposed nonlinear Poisson equation. In fact, physically, it is also very
diﬀerent as demonstarted in the solvation energy calculation. In general, to provide proper
prediction of experimental data, the nonlinear Poisson equation, as well as other smooth
dielectric Poisson equations, should have its own force ﬁeld parameters [153, 154, 37]. Note
that it is possible to set α = 1 so as to arrive at a parameter-free NLPE.

10

1.2.3

Nonpolar energy

The nonlinear Poisson equation provides a tool for electrostatic analysis. In general, the
electrostatic free energy of solvation is complemented by the nonpolar contribution. In the
present work, we adopt a popular nonpolar solvation model [114, 165, 37]

ρs U att (r)dr,

Gnp = γArea + pVol +

(1.19)

Ωs

where γ is the surface tension, Area represents the surface area of the molecule of interest,
p is the hydrodynamic pressure, Vol represents the volume occupied by the molecule, ρs is
the solvent bulk density, Ωs denotes the solvent accessible region and U att (r) is the solutesolvent van der Waals interaction potential at point r. The ﬁrst term γArea is the surface
energy, which describes the disruption of intermolecular and/or intramolecular bonds that
occurs when the surface of a molecule is created in the solvent. The second term pVol
measures the mechanical work of creating the vacuum of a biomolecular size in the solvent.
The hydrophobic eﬀect in the ﬁrst two terms are partially compensated by the third term
Ω s ρs U

att (r)dr,

which describes the attractive dispersion eﬀects near the solvent-solute

interface.
The nonpolar energy in Eq. (3.44) can be easily computed if the solvent-solute boundary
is explicitly given. However, in the present nonlinear Poisson model, there is no sharp interface between the solvent and the solute, which avoids the surface singularity and numerical
instability. To obtain a solvent-solute boundary, we deﬁne a linear mapping

S(ε) =

1
(εs − ε(| φ(r)|)) ,
εs − εm

11

(1.20)

where S(ε) is a solute characteristic function such that the solute domain Ωm is deﬁned
by S(ε) ≥ 0. Similarly, 1 − S(ε) is a solvent characteristic function such that the solvent
domain Ωs is deﬁned as 1 − S(ε) ≥ 0. Obviously, in the region εm < ε(| φ|) < εs , the
solvent domain Ωs overlaps the solute domain Ωm . Figure 1.2 illustrates the proﬁle of S(ε)
and 1 − S(ε) computed from Eq. (1.20). The impact of solvent polarity to the solvent cavity
or volume was discussed by Papazyan and Warshel [132].

1.2

S
1−S

1

0.8

0.6

0.4

0.2

0
−6

−4

−2
0
2
Y (angstrom)

4

6

Figure 1.2: The cross line proﬁles of S and 1 − S along the y-axis at x = 0 and z = 0 for
one-atom system.

In order to evaluate the nonpolar energy according to Eq. (3.44), one still needs to
compute the surface area and its enclosed volume. Once the surface proﬁle S is provided, it
is obvious that the volume of the solute can be expressed as

Vol =

dr =
Ωm

S(ε)dr,
Ω

12

(1.21)

where Ω is the total domain of the solvation problem. Additionally, we rewrite the surface
energy term as a function of the surface proﬁle S, which embeds a two-dimensional (2D)
surface in the R3 . By using the the coarea formula of the geometric measure theory [57, 180],
the surface area can also be expressed as a volume integral

| S(ε)|dr.

Area =

(1.22)

Ω

It is noted that Area only has contribution from

S = 0, where there is a transition region of

the solvent-solute boundary. Finally, the van der Waals dispersion term can be reformulated
as
ρs (1 − S(ε))U att (r)dr.

ρs U att (r)dr =
Ωs

(1.23)

Ω

We adopted the Weeks-Chandler-Andersen (WCA) decomposition based on the original WCA theory [176]. Therefore, the nonpolar solvation free energy Gnp can be deﬁned via the
surface proﬁle S(ε) given in Eq. (1.20),

γ| S(ε)| + pS(ε) + ρs (1 − S(ε))U att (r)dr.

Gnp =

(1.24)

Ω

For the attractive dispersion potential U att , which describes the solvent-solute coupling,
we have U att =

att
i Ui ,

where Uiatt accounts for the attractive van der Waals interaction

due to the ith solute atom. The Lennard-Jones (L-J) potential is one way to model the van
der Waals interaction. Here we consider the following 6-12 L-J potential at position r:

UiLJ (r)

= εi

σi + σs 12
−2
|r − ri |

13

σi + σs 6
,
|r − ri |

(1.25)

where εi characterizes the well-depth, σi and σs are ith solute atomic and solvent radii,
respectively, It is seen that UiLJ (r) achieves the minimum −εi at the position r = σi +
σs . There are diﬀerent ways to decompose L-J potential into attractive and repulsive
components. One is the ”6-12” decomposition, which is :

att,6/12
Ui
(r)
rep,6/12

Ui

= −2εi

(r) = εi

σi + σs 6
,
|r − ri |

(1.26)

σi + σs 12
.
|r − ri |

(1.27)

Another popular way is the Weeks-Chandler-Andersen(WCA) decomposition based on the
original WCA theory, which has the form of



 −εi ,
att,WCA
Ui
(r) =
 LJ
 U (r),
i

 LJ
 U (r) + εi ,
i
rep,WCA
Ui
(r) =

 0,

|r − ri | < σi + σs ;

(1.28)

|r − ri | ≥ σi + σs ,
|r − ri | < σi + σs ,

(1.29)

|r − ri | ≥ σi + σs .
att,WCA

In current work, we adopt the WCA decomposition, i.e. we set Uiatt (r) = Ui

(r).

It is worthwhile to point out that the distance in the present 6-12 Lennard Jones potential
in Eq. (1.25) are not exactly the same as those in the standard version due to the continuum
representation of solvent in our model. As a matter of fact, |r−ri | in Eq. (1.25) is the distance
between a speciﬁc position r in the solvent area and the center of a solute atom at ri , which
is no longer the distance between the center of a solute atom and that of a solvent atom in the
standard expression. Consequently, the setting of well depth εi diﬀers from those in popular

14

force ﬁelds such as from AMBER or OPLS . However the performance of the L-J potential
should be based on the same principle, i.e, the value of the L-J potential in the solvent caused
by a solute atom only depends on the distance from the center of the atom. This implies
that the value of L-J potential caused by a solute atom should be a constant on the van der
Waals surface of the atom. In other words, UiLJ (r) = εi

σi +σs 12
|r−ri |

σ +σ

i
− 2 |r−r s
i|

6

= Di if r

is on the vdW surface of the ith atom. Theoretically, the constant Di should have diﬀerent
values for diﬀerent types of atoms. For simplicity we use a uniform constant D to determine
the value of εi given σs and σi . In the present calculation, we set D = 1.0. Note that WCA
expression is chosen as the attractive van der Waals potential.
Thus, the above expression in Eq. (1.24) provides a practical way to compute the nonpolar solvation free energy based on the proposed NLPE (1.14).

1.2.4

Numerical Methods

1.2.4.1

Solution of nonlinear Poisson equation

Because the nonlinear Poisson equation (NLPE) admits a smooth dielectric proﬁle, it is
convenient to use the standard second order center diﬀerence scheme for the spatial discretization. Denote the position (xi , yj , zk ) by the pixel (i, j, k), then the NLPE can be
discretized as

1
1
ε(xi + h, yj , zk )[φ(i + 1, j, k) − φ(i, j, k)] + ε(xi − h, yj , zk )[φ(i − 1, j, k) − φ(i, j, k)]
2
2
1
1
+ε(xi , yj + h, zk )[φ(i, j + 1, k) − φ(i, j, k)] + ε(xi , yj − h, zk )[φ(i, j − 1, k) − φ(i, j, k)]
2
2
1
1
+ε(xi , yj , zk + h)[φ(i, j, k + 1) − φ(i, j, k)] + ε(xi , yj , zk − h)[φ(i, j, k − 1) − φ(i, j, k)]
2
2
= −4ρm (i, j, k)/h,

15

(1.30)

where h is the grid step size, φ(i, j, k) and ρm (i, j, k) are the discretized electrostatic potential
and molecular charge at grid point (xi , yj , zk ) respectively. In our work, the second order
interpolation (i.e., the trilinear mapping) is used to distribute charge ρm to the grid point
(i, j, k). The proﬁle of permittivity ε is obtained from the discretization as below

{ε(| φ|)}ijk = εm +

εs − εm
1 + 2kα T {| φ|2 }ijk
B

or
{ε(| φ|)}ijk = εm + (εs − εm )e

−α

n,

n = 1, 2,

{| φ|2 }ijk
2kB T ,

(1.31)

(1.32)

where
{| φ|2 }ijk = {φ2 }ijk + {φ2 }ijk + {φ2 }ijk ,
x
y
z

(1.33)

{φx }ijk = (φ(i+1)jk − φ(i−1)jk )/2h,

(1.34)

{φy }ijk = (φi(j+1)k − φi(j−1)k )/2h,

(1.35)

{φz }ijk = (φij(k+1) − φij(k−1) )/2h.

(1.36)

with

and

As a result, one attains a system of highly complex and nonlinear equation if | φ|2 for
ε is directly discretized. To reduce the computational cost, we solve the system by coupling
Eq. (1.30) and Eq. (1.31) or Eq. (1.32).
The far ﬁeld condition similar that used in the Poisson-Boltzmann equation [67] is employed for the nonlinear Poisson equation. As a result, the discretized nonlinear Poisson
16

equation can be cast into a linear algebraic equation system based on Eq. (1.30).

1.2.4.2

Initialization of dielectric proﬁle ε

To compute the dielectric proﬁle ε, we ﬁrst set up an initial value of ε as

ε(x, y, z) =



 εm , (x, y, z) ∈ D

 ε ,
s

(1.37)

otherwise

where the domain D is deﬁned as the region enclosed by the solvent accessible surface, i.e.,
D =

Nm
i=1 {r

: |r − ri | < ri + rp }. Here rp is the probe radius, ri is the atomic radius

and ri = (xi , yi , zi ), i = 1, · · · , Na are the atom centers of the ith atom and Nm is the
total number of the atoms for a given molecule. In the present work, ε(x, y, z) is ﬁxed as
εs outside the solvent accessible surface, and as εm within the van der Waals surface, i.e.,
(x, y, z) ∈

Na
i=1 {r

: ri < |r − ri | < (ri + rp )}, dielectric proﬁle ε(x, y, z) is determine by Eq.

(1.15) or Eq. (1.16). Although the initial value of ε has a sharp jump at the interface of the
solvent accessible surface, the iteration procedure of the solution of the nonlinear Poisson
equation will generate a smooth ε proﬁle connecting the solvent and the solute regions.

1.2.4.3

Evaluation of solvation free energy

In this subsection, we present the method for the evaluation of solvation free energies. Since
the solvation free energy is the energy released or absorbed when the solute molecule is
transferred from vacuum to the solvent environment, it is computed by

∆G = G − G0 ,

17

(1.38)

where G is the free energy calculated from the solvent environment while G0 is from the
vacuum. Typically , solvation free energy consists of polar and nonpolar contributions. Note
that the vacuum environment does not contribute to the nonpolar energy part. Therefore,
we have
∆G = Gp + Gnp − G0 = ∆Gp + Gnp ,

(1.39)

where ∆Gp = Gp − G0 and Gp and Gnp denote the polar and nonpolar contribution of G.
In the present work, the electrostatic solvation free energy ∆Gp is calculated by
N
1 m
Qi (φ(ri ) − φ0 (ri )),
∆Gp = Gp − G0 =
2

(1.40)

i=1

where φ and φ0 are electrostatic potentials in the solvent and the vacuum environment
respectively, and Qi is the partial charge on ith atom located at position ri . Finally, the
nonpolar solvation free energy is computed by Eq. (1.24).

1.2.4.4

Iteration procedure

The solution of the nonlinear Poisson equation requires appropriate iterations. The procedure
of coupling the Poisson equation and dielectric function is follows.
(Step 0): The workﬂow begins with the initialization of ε.
(Step 1): Then we solve NLPE for the electrostatic potential φ.
(Step 2): Next, the gradient

φ and thus ε is computed.

(Step 3): The electrostatic free energy ∆Gp is calculated.
(Step 4): We check the convergence by comparing |∆∆Gp |, which is the diﬀerence between
updated electrostatic free energy ∆Gnew and previous value ∆Gold , with a given tolerance.
p
p
If convergence is not reached, go back to (Step 1) until convergence. Otherwise, we stop the
18

iteration and compute ∆G, including both polar energy ∆Gp and s non-polar energy Gnp .
An overview of the algorithm is showed in the ﬂowchart in Figure 1.3.
Initialize ε

Solve NLPE for φ

Get

φ and ε

Obtain ∆Gp

New loop

no

|∆∆Gp | < τ ?

yes

Output ∆G

Figure 1.3: The ﬂowchart for the solution of the nonlinear Poisson equation.

1.3

Results and Discussion

The Poisson equation has been extensively applied in implicit solvent models to the solvation
analysis, which provides a computationally less expensive alternative to the integral equation
theory [18, 127], nonlocal dielectric theories [109] and density functional theory models [200].
In our work, we consider solvation applications of our nonlinear Poisson equation.

19

∆Gp (kcal/mol)

−50

−100

−150

−200
0

10

20

30

40

50

N

Figure 1.4: The convergence of electrostatic solvation free energies ∆Gp over number of
iteration (N) for one-atom system.

1.3.1

Validation by one atom system

We ﬁrst consider an atom of unit van der Waals radius and unit charge in the solvent. The
atomic center and the charge are both located at the origin. We set εm = 1 and εs = 80 for
dielectric constants in the solute and the solvent, respectively. A ﬁne mesh size of h = 0.25
˚ is used in our computation. The convergence of our computational algorithm with respect
A
to the number of iteration is shown in Figure 1.4 while the convergence over the grid size h
is validated as in Table 1.1.
Table 1.1: The convergence of ∆Gp with respect to grid size h for one-atom system.
grid size h
0.5
∆Gp
-183.8562

0.25
-183.6532

0.125
-183.6531

0.0625
-183.6531

In general, the dielectric function ε has a smooth proﬁle as shown in Figure 1.1. It
is interesting to note that as α → 0, the dielectric proﬁle converges to the discontinuous
dielectric constants of the classic linear Poisson model. It is worthwhile to further explore
20

20
Nonlinear Poisson model
Linear Poisson model

18

Volume (angstrom3)

16
14
12
10
8
6
4
2
0

0

2

4

α

6

8

10

Figure 1.5: The volume of the nonlinear Poisson model at a wide range of α values for
one-atom system with atomic radius 1 ˚, εm = 1 and εs = 80.
A
the impact of the hyperpolarization by using scaling factor α. To this end, we compute the
enclosed volume of the atom by Eq. (1.21) under a wide range of α values. In general,
the enclosed volume is larger than that of the linear Poisson model as shown in Figure 1.5.
Indeed, the volume converges to that of the linear Poisson model asymptotically α → 0.
In fact, the surface area computed with Eq. (1.22) converges to that of the linear Poisson
model as α → 0.
Having established the recovery of the linear Poisson model in terms of dielectric function,
surface area and enclosed volume, we now investigate polar solvation energy ∆Gp of the
nonlinear Poisson model. Figure 1.6 plots the behavior of ∆Gp for various α values. Unlike
other quantities, the solvation energy does not recover that of the linear Poisson model,
21

−100
Nonlinear Poisson model
Linear Poisson model

−110
−120

∆Gp (kcal/mol)

−130
−140
−150
−160
−170
−180
−190
−200

0

2

4

α

6

8

10

Figure 1.6: The electrostatic solvation energy of the nonlinear Poisson model at a wide range
of α values for one-atom system with atomic radius 1 ˚, εm = 1 and εs = 80.
A
i.e., the Kirkwood model [105], as α → 0. The plot shows that the Kirkwood model gives
∆Gp = −163.85 kcal/mol, while, the present nonlinear Poisson model predicts −183.65
kcal/mol. This behavior reminds us that the present nonlinear Poisson model is inherently a
smooth-permittivity model. It is well-known that without appropriate reparameterization,
smooth-permittivity or smooth dielectric models lead to an under-estimation of the polar
solvation energy [153, 154, 37]. It has been shown that the solvation energy of sharp dielectric
models can be recovered by reparameterizing atomic radii with a common factor [154, 37].
This in fact deﬁnes their own force ﬁeld parameters. In the present work, we perform a
similar parameter ﬁtting by adjusting the atomic radii so that the ∆Gp of nonlinear Poisson
model is the same as that of the Kirkwood model at the limit of α → 0. In this manner,
22

we can systematically determine the radii of commonly used atoms for our model. Table 1.2
lists the result of our radius set. We present both the enlarged radii as wells as the ratios for
the radius enlargement. Here the original radii are from the ZAP-9 force ﬁeld by Nicholls
etc [129]. It should be noted that values of atomic radius from diﬀerent models may diﬀer
signiﬁcantly. For example, the radius for Hydrogen is 1.20 ˚ from the Bondi’s compilation
A
but is adjusted to 1.10 ˚ by Nicholls [129]. Furthermore, since the ratios in Table 1.2 vary
A
from 1.06 to 1.12, it is more convenient to apply a common ratio to obtain our parameter set.
Figure 1.7 shows that a common ratio of 1.09 ﬁts well with all data. This is consistent with
our earlier ﬁnding in the diﬀerential geometry based solvation model [37] that a common
scaling factor of 1.1 gives rise to a good parameterization.
Table 1.2: Comparison of atomic radii and adjusted values (˚).
A
atom-type
H
N
O
Cl
O
S
F

1.3.2

radius [129]
1.10
1.40
1.76
1.82
1.87
2.15
2.40

adjusted-radius
1.23
1.54
1.88
1.95
2.00
2.30
2.54

ratio
1.12
1.10
1.07
1.07
1.07
1.07
1.06

Application to a set of 17 small molecules

To further validate the proposed nonlinear Poisson model, we apply it to the solvation analysis of a set of 17 small compounds. Various approaches, including quantum mechanics and
Poisson-Boltzmann theory, have been applied to this test set by Nicholls et al [129]. This
test set has also been employed in our earlier work to validate our diﬀerential geometry

23

Enlarged radius

3

2.5

2

1.5

1
1

1.5

2

2.5

Radius

Figure 1.7: The rescaled atomic radii for the one-atom system with ﬁxed ε = εm within the
solute domain.
based solvation models (DGSMs) [37, 38]. This test set is frequently used because its experimental data of solvation free energies are available. Furthermore, this set was found to be
quite challenging to theoretical models because of its hyperpolarization and the existence of
polyfunctional or interacting polar groups, which causes strong solvent-solute interactions.
To compare with experimental data of solvation, both polar and nonpolar solvation energies are required. We utilize the nonlinear Poisson equation (1.14) with the anisotropic
coeﬃcients given in Eqs. (1.15) and (1.16) for the electrostatic analysis. The nonpolar energy is computed according to Eq. (1.24), which in turn depends on the nonlinear Poisson
model for the evaluation of surface function S, volume, and surface area.
The structure and charge parameters of 17 compounds are taken from those of Nicholls
et al [129]. In particular, atomic charges are adopted from the OpenEye-AM1-BCC v1 parameters [94]. Atomic coordinates and radii are based on ZAP-9 parameters [129], in which
certain types of radii are revised from Bondi radii to improve the agreement with experimental free energy. In our computation, we carry out the smooth-dielectric reparameterization of
24

(a)

(b)

0

−2

−2
NLPE (kcal/mol)

2

0

NLPE (kcal/mol)

2

−4
−6

−4
−6

−8

−8

−10

−10

−12

−10

−5
Exptl (kcal/mol)

−12

0

−10

−5
DGSM (kcal/mol)

0

Figure 1.8: Correlation of total solvation free energies of 17 compounds between the nonlinear
dielectric model and diﬀerent data. (a) Nonlinear dielectric model and experimental data;
(b) Nonlinear dielectric model and DGSM.

25

atomic radii by using individual radii listed in Table 1.2. Figure 1.8(a) depicts the correlation
between the current results and experimental data. The root mean square error (RMSE) of
computation results is 1.75 kcal/mol and the average error is 1.40 kcal/mol when the radii in
Table 1.2 are used. In comparison, the linear Poisson model oﬀers a RMSE of 1.87 kcal/mol
by Nicholls et al [129]. Their explicit solvent approach, which is much more expensive, reduces the RMSE to 1.71 ± 0.05 kcal/mol [129]. Therefore, present nonlinear Poisson model
provides a relatively good prediction of solvation energies for this set of molecules.
It is interesting to further compare the performance of the present nonlinear Poisson model with the more complicated DGSM [37], which relies on the Laplace-Beltrami equation to
provide the surface function and dielectric proﬁle. The RMSE of the DGSM is 1.76 kcal/mol,
which is very close to the present result. Indeed, since both models admit smooth dielectric proﬁles and employ similar nonpolar energy functionals, there is very good correlation
between their solvation predictions as shown in Figure 1.8(b).
Table 1.3 gives a summary of the computation results by the present nonlinear model
with a comparison to the experimental data. From Table 1.3, it is seen that the error for this
set of 17 molecules varies from −2.41 to 3.86 kcal/mol. However, the major errors are from
the calculation of two benzamide molecules. The RMSE will drop from 1.75 kcal/mol to 1.25
kcal/mol without these benzamide compounds. The error with benzamides may be related
to the charges of the carbonyl oxygens and tertiary nitrogens in the OpenEye-AM1-BCC v1.
This issue will be further explored in our future work.
Similar to the linear Poisson model, the proposed nonlinear Poisson model is also capable
of providing surface electrostatic analysis. Figure 1.10 (a) and (b) depict surface electrostatic
potentials of four compounds obtained at a given isosurface value S = 0.5. This information
can be used in various practical applications.
26

Table 1.3: Comparison of free energies (kcal/mol) for 17 compounds between nonlinear model
and experimental data.
Compound
glycerol triacetate
benzyl bromide
benzyl chloride
m-bis(triﬂuoromethyl)benzene
n,n-dimethyl-p-methoxybenzamide
n,n-4-trimethylbenzamide
bis-2-chloroethyl ether
1,1-diacetoxyethane
1,1-diethoxyethane
1,4-dioxane
diethyl propanedioate
dimethoxymethane
ethylene glycol diacetate
1,2-diethoxyethane
diethyl sulﬁde
phenyl formate
imidazole

Gnp ∆Gp
∆G
2.63 -13.02 -10.40
1.61 -5.01 -3.41
1.57 -5.16 -3.59
2.44 -3.20 -0.76
2.23 -9.57 -7.34
2.11 -8.01 -5.90
1.66 -4.32 -2.66
1.89 -8.51 -6.62
1.74 -4.62 -2.88
1.16 -5.71 -4.55
2.13 -8.23 -6.10
1.17 -4.67 -3.49
1.84 -8.86 -7.02
1.77 -4.37 -2.60
1.52 -2.44 -0.92
1.55 -8.03 -6.49
0.96 -11.66 -10.70

Exptl Error
-8.84 -1.56
-2.38 -1.03
-1.93 -1.66
1.07
-1.83
-11.01 3.67
-9.76
3.86
-4.23
1.57
-4.97 -1.65
-3.28
0.40
-5.05
0.50
-6.00 -0.10
-2.93 -0.56
-6.34 -0.68
-3.54
0.94
-1.43
0.51
-4.08 -2.41
-9.81 -0.89

Finally, we test the solvation predictions of alternative nonlinear dielectric functions. As
shown in Table 1.4, three nonlinear Poisson models based on diﬀerent dielectric functions
provide essentially the same prediction of solvation free energies. Therefore, the proposed
nonlinear Poisson models are not sensitive to the functional form of the nonlinear dielectric
proﬁle.
Table 1.4: Errors of solvation free energies for a set of 17 small compounds obtained from 3
nonlinear models based on adjusted radius in Table 1.2.

RMSE (kcal/mol)
Average error (kcal/mol)

n = 1 n = 2 exponential
model model
model
1.7505 1.7565
1.7531
1.4011 1.4051
1.4035

27

0
−500

MIBPB−III vs NLPE
DGSM vs NLPE

NLPE (kcal/mol)

−1000
−1500
−2000
−2500
−3000
−3500

−3000
−2000
−1000
MIBPB−III / DGSM (kcal/mol)

0

Figure 1.9: Correlation between the results obtained by the MIBPB-III [67] and the nonlinear
dielectric model for electrostatic solvation free energies of 20 proteins.

1.3.3

Application to a set of 20 proteins

After demonstrating the success in the validation and application in small compounds, we
examine our nonlinear Poisson model for the electrostatic solvation of biomolecules. Here
we chose 20 proteins which are frequently used in previous studies [67, 37]. The number of
atoms for this set of proteins ranges from 519 to 2809.
The initial structure data of all proteins are taken from the Protein Data Bank (PDB)
at http://www.pdb.org/. Based on the initial structure data, the hydrogen atoms which are
typically missing from the X-ray data are added to the structures and extra water molecules
that are attached to proteins are excluded so as to obtain full all-atom models. Partial charges

28

(a)

(b)

(c)

(d)

Figure 1.10: Surface electrostatic potentials of four molecules at their isosurface value
S = 0.5. (a) small compound: diethyl propanedioate; (b)small compound: n,n-dimethyl-pmethoxybenzamide; (c)protein: 1vjw; (d) protein: 2pde.
at atomic sites as well as atomic radii in angstroms are assigned from the CHARMM27 force
ﬁeld [120]. Note that a common radius scaling factor of 1.1 is used in this computation.
n-polar energy are the same as before.
The results for 20 proteins are summarized in Table 1.5 in the Supporting Information.
The magnitude of electrostatic solvation free energies for these proteins are much larger
29

0
−1
−2

NLPE (kcal/mol)

−3
−4
−5
−6
−7
−8
−9
−10
−11

−10

−8

−6
−4
Exptl (kcal/mol)

−2

0

Figure 1.11: Correlation of solvation free energies predicted by the NLPE model and experimental data for 21 compounds.
than that of the previous set of 17 small compounds due to multiple charges. Since the
experimental data for the total solvation free energy for this set of proteins is not available
yet, we compare the electrostatic solvation free energies with those from the DGSM [37] and
from the MIBPB-III [67]. It is seen that the results are close to each other. An illustration
of the correlation between the present nonlinear model and the DGSM is depicted in Figure
1.9. The results of three models are quite close. Figure 1.10 (c) and (d) display the proﬁle of
electrostatic potential on the isosurface value S = 0.50 for two proteins, i.e., 1vjw and 2pde,
which suggests that the present nonlinear Poisson model is also robust for the visualization
of electrostatic surface potentials.

30

Table 1.5: Electrostatic solvation free energies (kcal/mol) for 20 proteins by diﬀerent models.

PDB-ID
1ajj
2erl
1bbl
1vii
2pde
1sh1
1fca
1uxc
1fxd
1vjw
1bor
1hpt
1mbg
1r69
1neq
451c
1a2s
1svr
1a63
1a7m

1.3.4

No. of atoms
519
573
576
596
667
702
729
809
824
828
832
858
903
997
1187
1216
1272
1435
2065
2809

Gp
MIBPB-III
-1137.2
-948.8
-986.8
-901.5
-820.9
-753.3
-1200.1
-1138.7
-3299.8
-1237.9
-853.7
-811.6
-1346.1
-1089.5
-1730.1
-1024.6
-1913.5
-1711.2
-2373.5
-2155.5

(kcal/mol)
DGSM NLPE
-1178.5 -1109.8
-935.8 -930.7
-965.9 -956.8
-892.0 -878.3
-843.0 -850.9
-771.8 -797.6
-1200.6 -1223.1
-1125.7 -1109.4
-3291.9 -3307.0
-1226.6 -1233.3
-871.4 -872.3
-808.2 -793.0
-1328.2 -1317.9
-1072.7 -1055.1
-1713.9 -1694.1
-1020.6 -1019.1
-1900.3 -1910.7
-1704.6 -1682.6
-2380.5 -2384.5
-2179.8 -2201.0

Application to temperature eﬀects

It is well known that solvation free energies are temperature dependent [163]. Although
most implicit solvation models are developed for room temperature, temperature eﬀects
were considered in the SM6T model [30] and by Elcock and McCammon [56]. In this work,
we demonstrate the ability of the NLPE for modeling the temperature eﬀect on solvation
free energies. First, one notes that the NLPE explicitly depends on the temperature as
shown in Eq. (1.15). Additionally, dielectric constants, εm and εs , are known to vary with
temperature [115]. The dielectric constant of water reduces from 80.10 at 298K to 55.58 at
373K, and its values at other temperatures can be well approximated by a cubic polynomial

31

[115]. In implicit solvation models, εm is known to increase as the temperature increases
because of atomic ﬂuctuation [164]. In the present model, we set εm = 1 at 298K and εm = 2
at 373K, with a linear interpolation for other temperatures. It should also be noted that
one may need to carefully investigate temperature eﬀects on other parameters, namely, the
atomic radius and charge, in the NLPE to fully explore the present model. For simplicity,
atomic radius and charge are treated as temperature independent since their change with
respect to temperatures is neglectable [30].
To validate the NLPE model for temperature eﬀects, we consider a set of compounds
whose free energies of solvation at diﬀerent temperatures are available from experimental
measurement [30]. The structural data of these compounds can be computed by using PubChem (http://pubchem.ncbi.nlm.nih.gov/). However, the corresponding electronic charge
densities are obtained from a quantum approach developed by Chen and Wei [39].
Two diﬀerent temperatures are considered to illustrate temperature eﬀects on solvation
free energies of 21 compounds. Due to the availability of experimental data, considered temperatures are not the same for all the compounds, as shown in Table 1.7. As the temperature
increases, the average increase in solvation energies is 1.13 kcal/mol in experimental data
while that predicted by the NLPE model is 1.05 kcal/mol. The RMSE of solvation free energies predicted by the NLPE model is 1.23 kcal/mol at low temperatures and 1.20 kcal/mol
at high temperatures. Figure 1.11 illustrates the correlation between predicted solvation free
energies and experimental data at all temperatures. Obviously, the proposed NLPE model
does a good job in modeling the temperature eﬀect on the solvation free energy.
Note that by adjusting the dielectric constants, the classical Poisson equation can also be
used to model the temperature eﬀect. Table 1.6 shows the RMSE and average ∆∆G of two
models at low and high temperatures for 21 compounds. It is seen that the present NLPE
32

Table 1.6: Comparison of RMSE and ∆∆G predicted by the classical Poisson model and the
NLPE model for 21 compounds (The average ∆∆G of experimental data is 1.13 kcal/mol).
Model
RMSE (kcal/mol) at low T
RMSE (kcal/mol) at high T
Average ∆∆G (kcal/mol)

classical Poisson
1.56
1.65
0.86

NLPE (n=1)
1.23
1.20
1.05

model performs much better than the classical Poisson equation.

1.4

Chapter conclusion remarks

The classic Poisson theory neglects hyperpolarizations, which are important under a strong
electrical ﬁeld, or involving highly charged nonlinear inhomogeneous media. We propose a
new electrostatic solvation free energy functional to partially account for the eﬀect of hyperpolarizations to solvation analysis. A nonlinear Poisson equation is derived from the
Euler-Lagrange equation. The present nonlinear Poisson equation gives rise to a smooth dielectric function, which by-passes both conceptual and computational diﬃculty of employing
a sharp solvent-solute boundary. By using geometric measure theory, we also formulate a
nonpolar solvation model based on the dielectric proﬁle obtained from the nonlinear Poisson
model. The proposed solvation models are extensively validated with the Kirkwood model
and experimental measurements of 17 molecules. Our new solvation models out-perform
the classic Poisson equation based solvation model. In fact, its performance is very close to
the explicit solvation model, which is much more computationally expensive. Application of
the proposed nonlinear Poisson model is considered to electrostatic analysis of 20 proteins.
Additionly, the test results for the set of 21 compounds at diﬀerent temperatures further
validate our model.
33

Table 1.7: Temperature eﬀect for 21 compounds between the NLPE model and experimental
results.

Compd
anthracene
phenanthrene
1-methylphenanthrene
ﬂuoranthene
Acetaldehyde
1-Methoxypropane
tetrahydrofuran
2-Methoxypropane
ethyl-formate
formic-acid
1,2-ethanediol
naphthalene
1-methylnaphthalene
acenaphthylene
acetone
2-Butanone
cyclopentanone
benzylal
cyclohep
2methylc
cyclohexylm
Mean
RMSE

∆G at low T
T Exptl NLPE
284 -4.14 -4.75
277 -4.29 -5.67
277 -4.01 -5.51
283 -5.12 -5.04
273 -4.03 -3.51
273 -2.35 -2.15
293 -3.58 -2.31
283 -2.32 -2.16
278 -2.93 -2.86
278 -7.09 -5.68
298 -9.35 -8.19
277 -3.01 -1.13
277 -3.08 -5.06
277 -3.78 -2.71
274 -4.38 -4.64
283 -4.55 -3.61
273 -5.83 -7.31
273 -7.61 -5.04
293 -5.77 -5.35
273 -5.47 -6.73
303 -4.09 -5.30
1.23

34

∆G at high T
T Exptl NLPE
304 -3.45 -3.23
309 -3.54 -4.85
304 -3.58 -5.08
328 -3.86 -3.33
313 -3.06 -3.37
298 -1.64 -1.58
343 -2.72 -1.43
298 -1.90 -1.61
349 -2.25 -1.80
309 -6.74 -5.06
348 -8.39 -7.42
309 -2.13 -0.62
313 -2.06 -4.45
304 -2.93 -2.15
333 -3.17 -3.58
373 -2.46 -1.68
364 -4.31 -5.92
328 -5.89 -3.70
364 -3.59 -3.86
353 -3.04 -4.46
363 -2.47 -4.03
1.20

∆∆G
Exptl NLPE
0.69
1.52
0.75
0.82
0.43
0.43
1.26
1.71
0.97
0.14
0.71
0.57
0.86
0.88
0.42
0.55
0.68
1.06
0.35
0.62
0.96
0.77
0.88
0.50
1.02
0.61
0.85
0.56
1.21
1.07
2.09
1.93
1.52
1.40
1.72
1.34
2.18
1.48
2.43
2.27
1.62
1.27
1.13
1.05

Chapter 2
PDE transform for hyperbolic
conservation laws

2.1

Introduction

Hyperbolic systems of nonlinear conservation laws

ut + f(u)x = 0

(2.1)

u(x, 0) = u0 (x)

(2.2)

with an initial condition

have attracted great attention in the past few decades in mathematical, scientiﬁc and engineering communities due to their practical applications in ﬂuid mechanics, aerodynamics,
and nano-bio systems, to mention only a few. The solution to this class of problems may
not exist in the classical sense because of possible discontinuities in the initial condition,
material interface, singularity formation, turbulence, blow-up, etc.
Both global and local methods have been developed for hyperbolic conservation laws.
Many up-to-date local methods have been proposed for shock-capturing, turbulence and

35

shock interaction, including weighted essentially non-oscillatory (WENO) scheme [95], central schemes [20, 111, 116, 128], arbitrary-order non-oscillatory advection scheme [158], gas
kinetic [187], anisotropic diﬀusion [131], conjugate ﬁlters [86] and image processing based
algorithms [79, 179]. An important factor that contributes to the success of the above mentioned local schemes in the shock-capturing is their appropriate amount of intrinsic numerical
dissipation, which is introduced either by explicit artiﬁcial viscosity, upwinding, relaxation,
or by local average strategy in non-oscillatory central schemes [98]. Indeed, the characteristic
decomposition based on Roe’s mean matrix can also be considered as a local average of the
Jacobian matrix. The relation between some approximate Riemann solvers and relaxation
schemes was analyzed by LeVeque [113]. Local characteristic decomposition is not necessary in low-order methods because of intrinsic numerical dissipation, while it seems to be
indispensable in high-order methods [135]. In general, local and low order methods perform
well for problems whose Fourier responses of the solution focus predominantly in the low
frequency region. For this class of problems, ﬁrst order or second order Godunov type of
schemes can be very eﬃcient in balancing accuracy and eﬃciency. When local and low order
methods are used for resolving shocks in ﬂows with ﬁne structural details or highly oscillatory
patterns, their numerical dissipation is usually too large to oﬀer informative results.
Spectral methods, or global methods, on the contrary, produce little numerical dissipation and dispersion in principle when applied to approximate spatial derivatives. It is well
known that spectral methods are some of the most accurate and eﬃcient approaches for
solving partial diﬀerential equations (PDEs) arising from scientiﬁc and engineering applications [28, 62, 159]. Therefore, it is highly desirable to use spectral methods for the solution of hyperbolic conservation laws, because the multiscale features, including Kolmogorov
microscales, in the hyperbolic conservation law systems require high resolution methods.
36

Nevertheless, when spectral methods are applied to the approximation of spatial derivatives
of a discontinuous function, which often occurs in hyperbolic conservation law problems,
one encounters Gibbs’ oscillations [77]. Most previous investigations are aimed at improving
the rate of convergence away from the discontinuity while recovering smooth solutions from
the contamination of Gibbs’ oscillations. The suppression of Gibbs’ oscillations is necessary
in order to avoid unphysical blow-ups in the time integration. Therefore, it has been of
tremendous interest in modifying spectral methods for hyperbolic conservation law systems
in the past two decades [112, 151]. There are two general types of approaches in spectral
based methods for hyperbolic conservation law systems: (1) explicit artiﬁcial viscosity, e.g.
spectral viscosity method proposed by Tadmor [155], and (2) ﬁltering. It is expected that
the appropriate use of spectral methods enables us not merely to capture the shock, but
also to resolve the delicate features, immersed interface, and underlying ﬁne structures of
the ﬂow.
Filters are designed to apply either in the spectral domain, called spectral ﬁlters, or in the
physical domain. Typical spectral ﬁlters include Lanzos ﬁlter, raised cosine ﬁlter, sharpened
raised cosine ﬁlter, Krasny ﬁlter [110] and exponential cutoﬀ ﬁlter, as listed by Hussaini
et al. [93]. More sophisticated and eﬀective ﬁlters in spectral domain are Vandeven’s pth
order ﬁlter [162], and Gottlieb and Tadmor’s regularized Dirichlet function [78]. A ﬁlterbased Reynolds-averaged Navier-Stokes approach was developed to improve the predictive
capability considerably in comparison to the standard k − model [101] .
Filters in the physical domain are also commonly used alternatives to spectral ﬁlters. In
the framework of spectral methods, it is generally more diﬃcult to design appropriate ﬁlters
in the physical domain than in the Fourier domain. A simple procedure is to make use of
numerical dissipation contained in some high-order shock-capturing schemes [27], such as
37

the ENO scheme, where, actually, numerical dissipation was introduced both in the Fourier
domain (via an exponential ﬁlter) and in the physical domain (via ENO polynomial ﬁlter).
Such a strategy was employed by Yee et al. [190] to construct characteristic ﬁlters in the
framework of ﬁnite diﬀerence methods. Gegenbauer polynomials are used to resolve the
oscillatory partial Fourier summation [77]. Promising numerical results were generated by
using ﬁlter approaches [51, 52, 75].
About a decade ago, we proposed a conjugate ﬁlter oscillation reduction (CFOR) scheme
[85, 181, 205, 206, 151] for hyperbolic conservation laws. This scheme was constructed within
the framework of a local spectral wavelet method, namely, the discrete singular convolution
(DSC) algorithm [167, 177, 183]. Here, ‘conjugate ﬁlters’ means that the eﬀective wavenumber range of the low-pass ﬁlter is largely overlapped with that of the high-pass ﬁlter used
for the approximation of spatial derivatives. In fact, the DSC algorithm is used to behave
as both low-pass ﬁlters and high-pass ﬁlters. Extensive validation of the CFOR scheme
over a wide range of shock-capturing problems has been carried out [85, 181, 205, 206]. We
demonstrated that CFOR scheme provides some of the highest grid resolution, i.e., 5 points
per wavelength (PPW), for the interaction of shock and entropy waves, and for many other
challenging problems involving natural high frequency oscillations [205, 206, 151].
Despite great eﬀort in the past few decades, the eﬃcient application of ﬁlters for hyperbolic conservation law systems remains a challenging problem. To design eﬃcient ﬁlter
methods, one must control a number of ﬁlter properties, such as ﬂatness, ripple, ﬁlter length,
eﬀective frequency range and length of transition band, to name only a few. Normally, it
is desirable to utilize ﬁlters that are free of dispersion errors, ﬂat while having very small transition band, short in length while having high resolution. Obviously, some of these
properties are conﬂicting with each other. Adjustable parameters have to be employed to
38

tune ﬁlter design properties in the application. In addition to the diﬃculties in controlling
ﬁlter properties, there are intrinsic mathematical challenges in developing ﬁltered spectral
methods. First, the solutions to diﬀerent hyperbolic conservation law equations may have
diﬀerent Fourier spectral distributions. Additionally, even for a single conservation law equation, the characteristic of the Fourier spectral distribution may change during the time
integration. Finally, the same conservation law equation can exhibit dramatically diﬀerent
Fourier spectral distributions when it is applied to diﬀerent physical problems, i.e., problems
with diﬀerent initial, interface and boundary conditions. Therefore, an optimal ﬁlter has to
be fully adaptive to hyperbolic equations, to initial, interface and boundary conditions, and
to the variation of spectral characteristic during the time integration. Unfortunately, such a
ﬁlter approach still does not exist yet. Given the complexity in the hyperbolic conservation
law and in ﬁlter properties, it is unlikely that there will be a perfect solution to these challenges in the near future. Consequently, these challenging and important problems call for
the further study of spectral ﬁlter approaches.
To put the matter in perspective, currently, there is no such a method that is perfect
for all of the tasks in hyperbolic conservation law systems. The challenges are due to the
intrinsic mathematical diﬃculties of hyperbolic conservation law equations and the inherent
limitation in the understanding of the fundamental physics involved in the hyperbolic conservation laws as noticed by Nobel Laureate Richard Feynman. As it is well known, hyperbolic
conservation laws can be formally derived from the (quantum) Boltzmann equation, which
in turn can be derived from the (quantum) Liouville equation equation via the BBGKY hierarchy [179]. Many physical assumptions are built into these derivations and their validities
may not be warranted under certain physical situations. There is a pressing need to revisit
basic assumptions and fundamental physics for the hyperbolic conservation laws for complex
39

physical applications.
Most recently, we have introduced PDE transform as a new approach for the analysis
of signals, images, surfaces and data [171, 169, 203]. The PDE transform is based on a
family of arbitrarily high order nonlinear PDEs ﬁrst introduced by Wei for edge-preserving
image restoration in 1999 [178] and PDE based high-pass ﬁlters proposed by Wei and Jia
in 2002 [182]. Variational models of PDE transform have also been proposed [169]. With
an iterative procedure [170, 171] to incorporate appropriate residues, i.e., initial conditions,
the PDE transform is able to extract functional mode functions (FMFs) and allow perfect
reconstruction. By FMFs, we mean the mode components which share same band of frequency distribution as well as same category of physical functions, i.e., trend, edge, texture,
feature, trait, noise etc. Using the FMFs obtained from the PDE transform, secondary
processing, or post-processing, can be performed to achieve desirable tasks, such as edge
detection, trend estimation, image enhancement, denoising, texture quantiﬁcation, segmentation, feature extraction, pattern recognition, etc. The PDE transform can perform as
tunable ﬁlters. By adjusting the order of the PDE transform, i.e., the highest order of the
PDE, one controls time-frequency localization, while by adjusting the diﬀusion coeﬃcients
or propagation time, one obtains desirable frequency precision or multiresolution analysis
[203]. The PDE transform can be tuned according to the input data to provide desirable
mode information. The full process of the PDE transform is nonlinear even if linear PDEs
are employed. Unlike the Fourier transform or the wavelet transform, the PDE transform
conserves the data representation. The resulting functional modes are still in the original
data form. The PDE transform has been applied to image analysis [171, 169], signal processing [171] and biomolecular surface construction [203]. High order PDE transforms with
their order being much higher than 4 are found to play a vital role in signal analysis, image
40

processing and surface generation.
The objective of present work is to explore the utility of the PDE transform for the
solution of hyperbolic conservation laws. In fact, such a utility can be realized in a number
of ways. For simplicity, we make use of a fast PDE transform, i.e., the PDE transform realized
by the fast Fourier transform in a single time stepping. Such a fast algorithm is paired with
the Fourier pseudospectral method (FPM) for solving hyperbolic conservation law equations
and suppressing Gibbs’ oscillations. Therefore, the fast PDE transform implemented in the
present work is essentially a spectral ﬁlter approach. Nevertheless, the PDE transform can
be implemented with geometric structures using mean-curvature ﬂow and/or Willmore ﬂow
type of nonlinear PDEs. However, this aspect is beyond the scope of the present work. In
fact, the PDE transform has an adjustable eﬀective wavenumber range that makes it viable to
capture ﬁne ﬂow structures, which is a desirable objective of spectral methods for hyperbolic
conservation law systems. It is this adjustable eﬀective wavenumber range that controls
the resolution of the overall scheme. In our design, this adjustable eﬀective wavenumber
range can be varied by the highest order of the PDE transform and the duration of time
propagation according to the resolution requirement of a hyperbolic conservation equation
and the physical problem of interest. Diﬀerent orders of the PDE transform have diﬀerent
magnitude responses and adjustability in the Fourier domain, which in turn inﬂuences the
accuracy and resolution of the PDE transform based FPM. It is this ﬂexibility that makes
the present method applicable to a wide variety of hyperbolic conservation law systems. The
performance of the proposed method is extensively validated and compared with those of
other approaches in the literature.

41

2.2

Theory and algorithm

To establish notation and enhance the basic understanding of the proposed PDE transform
strategy for systems of hyperbolic conservation laws, we present a brief introduction of arbitrarily high order nonlinear PDEs and PDE transforms. The detailed numerical algorithm
for time integration of evolution equations of hyperbolic conservation laws is described.

2.2.1

Arbitrarily high order nonlinear PDEs

Last two decades have witnessed a dramatical growth in PDE based methods for image
analysis. Using the diﬀusion equation for image denoising was pioneered by Witkin in 1984
[186]. Although, Witkin’s algorithm is mathematically simple, it provides the foundation
for much later development. Since nonlinear PDEs and variational calculus are attractive
research topics, anisotropic diﬀusion equation introduced by Perona and Malik [134] and
total variation models pioneered by Rudin, Osher, and Fatemi [141] have had much impact
in the applied mathematical community. These earlier approaches are based on the use of
the second order PDEs for image or surface analysis. In 1999, Wei introduced the ﬁrst family
of arbitrarily high order nonlinear PDEs for edge-preserving image restoration [178]

∂u(r, τ )
=
∂τ
q
where r ∈ Rn ,

· dq (u(r, τ ), | u(r, τ )|, τ )

2q u(r, τ )

+ e(u(r, τ ), | u(r, τ )|, τ ), (2.3)

∂
= ∂r , u(r, τ ) is the processed image function, dq (u(r, τ ), | u(r, τ )|, τ )

are edge sensitive diﬀusion coeﬃcients and e(u(r, τ ), | u(r, τ )|, τ ) is enhancement operator.
Equation (2.3) is subject to the initial image data u(r, 0) = X(r) and appropriate boundary
condition. Its construction was motivated by the super ﬂux in the pattern formation in

42

alloys, glasses, polymer, combustion, and biological systems and the Perona-Malik equation
[134]. Its essential idea is to accelerate the noise removal by higher order derivatives, which is
more eﬃcient in noise dissipation. An interesting feature is that one can recover the PeronaMalik equation by setting q = 0 and e(u(r, τ ), | u(r, τ )|, τ ) = 0. Therefore, Eq. (2.3) is also
called generalized Perona-Malik equation. The fourth order version of Eq. (2.3) has been
applied to image denoising and restoration by many researchers [178, 119, 70, 69].
As in the original Perona-Malik equation, one can choose the hyperdiﬀusion coeﬃcients
dq (u, | u|, τ ) in Eq. (2.3) in many diﬀerent ways. One of the popular choices is the Gaussian
form

dq (u(r, τ ), | u(r, τ )|, τ ) = dq0 exp −

| u|2
,
2
2σq

(2.4)

where the values of constant dq0 depend on the noise level, and σ0 and σ1 were chosen as
the local statistical variance of u and

u

2
σq (r) = | q u −

q u|2

(q = 0, 1).

(2.5)

The notation Y (r) above denotes the local average of Y (r) centered at position r. In this
algorithm, the measure based on the local statistical variance is important for discriminating image features from noise. As such, one can bypass the image preprocessing, i.e., the
convolution of the noise image with a test function or smooth mask in the application of the
PDE operator.
In the past decade, high order nonlinear PDEs, particularly fourth order nonlinear PDEs,
have attracted much attention in image analysis [31, 178, 32, 191, 157, 81, 119, 82, 13]. Unlike

43

the classical second order PDEs, higher order PDEs are able to suppress high frequency
oscillations, including noise, at much faster rates. Mathematical analysis of high order
nonlinear PDEs has been a popular research topic in applied mathematics. Bertozzi and
Greer proved the existence and uniqueness of the solution of fourth order edge-preserving
PDEs in Sobolev space with H 1 initial data and a regularized operator [19, 81, 82]. A similar
mathematical analysis was carried out by Xu and Zhou [188] for fourth order nonlinear PDEs.
Recently, Jin and Yang have presented an interesting result in which they show that the
mathematical structure of Wei’s fourth order nonlinear PDE diﬀers from other fourth order
PDEs derived from the variational formulation [99]. They have proved the existence of the
strong solution of such a PDE [99]. Because high order PDEs are subject to strict stability
constraints in their numerical solutions, the development of stable and eﬃcient numerical
techniques for higher order PDEs is an important issue, except for digital image processing
where the grid size is normally unit. Witelski and Bowen proposed alternating direction
implicit (ADI) schemes to solve arbitrarily high order nonlinear PDEs [185]. Similar ADI
schemes were also constructed in our work [13].
Recently, another family of arbitrarily high order geometric PDEs was proposed for
surface formation and evolution with application to biomolecular systems [13],

∂S
= (−1)q
∂τ

g(|

2q S|)

where S is the hypersurface function, g(|

( 2q S)

·

g(|
2q S|)

2q S|)

=1+|

+ P (S, | S|),

2q S|2

(2.6)

is the generalized Gram

determinant and P is a generalized potential term, including microscopic potential eﬀect in
biomolecular surface construction. Equation (2.6) was designed as a generalization of other
important geometric PDEs. For example, when q = 0 and P = 0, it reduces to the mean

44

curvature ﬂow equation used in our earlier construction of minimal molecular surfaces [16],
while when q = 1 and P = 0, it is a surface diﬀusion ﬂow [13]. Certainly, Eq. (2.6) can also
be regarded as a variant of Wei’s earlier arbitrarily high order PDE (2.3). It is interesting
to note that the molecular surface generated by Eq. (2.6) has a morphology distinguished
from that generated by using the second order geometric PDEs [13].
Variational approaches provide important theoretical formulations for physical and biological systems [16, 180, 37, 38]. Variational derivation of molecular surfaces are formulated
in many of our recent work [180, 37, 38, 34]. Similar approaches are commonly used in
image analysis in applied mathematics [140, 16, 14, 49]. Didas et al. discussed the variation
formulation of high order nonlinear PDEs [49]. Here we provide alternative expressions for
the PDE transform by variation. We construct the energy functional as
 
E(u,

u, 2 u, · · · , m u) =



m



| q u|2  + (X − u)2  dr,

Λ 

(2.7)

q=1

where

is a constant, X is the original data, (X − u)2 is the ﬁdelity term and Λ(·) is an

appropriate penalty function. The most commonly used penalty function is the Tikhonov
form

Λ(x2 ) = x2 .

(2.8)

Other useful forms include [49, 119]
1

Λ(x2 ) =

σ 2 + x2 2

2
2
Λ(x2 ) = ex /2σ .

45

(2.9)
(2.10)

Energy functional (2.7) can be minimized by using the Euler-Lagrange equation,


m
q

q

q

| q u|2 

Λuq 

q=1

where



m

qu +

(X − u) = 0

(2.11)

q=1

is a normal product for even q and an inner product for odd q. Here

Λuq = (−1)q+1 ∂Λ/∂| q u|2 .

(2.12)

By solving Eq. (2.11), one realizes the energy minimization. A practical way to solve Eq.
(2.11) is to introduce an artiﬁcial time τ and convert Eq. (2.11) into a time dependent PDE


m
q

∂τ u =
q=1

q



m

| q u|2 

Λuq 

qu +

(X − u).

(2.13)

q=1

Equation (2.13) is essentially equivalent to the earlier PDE transform [169].

2.2.2

PDE transform

One of the important properties of the PDE transform is its ability to extract mode functions
from a given data X. To illustrate this point, we denote the solution of Eq. (2.3) and/or
ˇ
Eq. (2.13) as X(τ ) such that

ˇ
X k (r, τ ) = LX k (r)

46

(2.14)

ˇ
where L is a low-pass PDE transform satisfying Lu(r, 0) = u(r, τ ), X k (r, τ ) are mode functions and X k (r) is the kth residue function deﬁned by

X 1 = X(r)

(2.15)

and
k−1

Xk

=

X1

ˇ
Xj,

−

∀k = 2, 3, · · · .

(2.16)

j=1

Obviously, there is a perfect reconstruction of the original data X in terms of all the mode
functions and the last residue
k−1

ˇ
Xj + Xk.

X=

(2.17)

j=1

Note that the PDE transform given in Eq. (2.14) recursively extracts mode functions based
on the input residue function. This procedure is nonlinear even if a linear PDE operator is
used because the initial value changes during the repeated operations.
The ﬁrst mode produced by the PDE transform described above is the trend of the
data. The residue of the trend is an edge function, including possible noisy components.
By systematically repeating the low-pass PDE transform (2.14), one can extract all the
desirable higher order mode functions. In our earlier work, high-pass PDE transforms were
also constructed in which the ﬁrst mode is edge type of information or possible noise; while
the ﬁnal residue is the trend [171].
Solving arbitrarily high order PDEs, such as Eq. (2.3) and/or Eq. (2.13) in the PDE
transform can be a very diﬃcult issue for some practical application. A main diﬃculty is
47

the stability constraint as the time stepping is normally proportional to the 2mth power of
spatial grid spacing, where 2m is the highest order of the PDE transform. An exception is
digital image processing, in which the spatial grid spacing is usually unit and thus bypasses
the stability constraint. However, spatial grid spacing is normally smaller than one in most
other applications.
There is another way that one can entirely bypass the stability constraint and obtain the
solution in a single time step. This approach relies on the use of the fast Fourier transform
(FFT) for the solution of evolution PDEs [151]. To this end, we linearize Eqs. (2.3) and
(2.13), and assume the form
m

∂τ v =

(−1)j+1 dj 2j v + (X k − v),

τ ≥ 0,

(2.18)

j=1

where dj > 0, ∼ 0 and X k ∈ Rn is the kth residue of the data. When Eq. (2.18) is subject
to initial value v(r, 0) = X k and periodic boundary conditions, it can be easily solved with
the Fourier pseudospectral method. The solution algorithm developed in the above Fourier
domain is called a fast PDE transform. In the present work, we explore the use of the fast
PDE transform for integrating hyperbolic conservation law systems.

2.2.3

Numerical algorithm

Both linear and nonlinear PDE transforms can be used for systems of hyperbolic conservation
laws. For simplicity, we consider only linear PDEs in the present work. The FPM is used
for the basic discretization of hyperbolic conservation law equations, while the fast PDE
transform is implemented as spectral ﬁlters.

48

The numerical scheme can be expressed as the following two-step procedure:

¯
U k+1 = BU k

(2.19)

¯
U k+1 = LU k+1 ,

(2.20)

where the operator B in Eq. (2.19) is the basic time integration from time k to time k + 1
and the operator L in Eq. (2.20) is a possible application of the PDE transform. Here, U k
and U k+1 is the numerical solution of u in Eq. (2.1) at time step k and k + 1, respectively,
¯
while U k+1 is the intermediate numerical solution of u after implementing operator B on
U k.
The operator B in Eq. (2.19) can be deﬁned by a general numerical method for the time
evolution. Here, we use the fourth-order Runge-Kutta scheme. The FPM is utilized for the
spatial discretization of f (u)x in the frequency domain.
Speciﬁcally, after carrying out the fast Fourier transform on f , we obtain the frequency
response f , i.e. f = F F T (f ). Then iω f (ω) is the frequency response of fx , where i =

√

−1

is the imaginary root. By doing the inverse fast Fourier transform on iω f (ω), we obtain the
spatial discretization of f (u)x as IF F T (iω f (ω)).
For simple and continuous problems, scheme (2.19) works well. However, for hyperbolic
conservation systems involving discontinuity, the accumulation of Gibbs’ oscillations as time
evolves may result in spurious solutions or even numerical blow ups. Thus, we make use of
PDE transform [151] to eliminate the possible Gibbs’ oscillations from basic time integration
(2.19). The application of the PDE transform is controlled (i.e., turned on or turned oﬀ) by
an adaptive sensor. The sensor is characterized by a measure of high frequency denoted by
µ

and regulated by a threshold value η. Once the increment in high frequency measure,

49

∆µ, exceeds the threshold η, the PDE transform in Eq. (2.20) is implemented. In our test,
we apply a straightforward high frequency measure

M

as the TVD sensor, which is

deﬁned by

¯
µ(U k+1 )

¯ k+1 ¯
| Ui+1 − Uik+1 |,

=

(2.21)

i

¯
¯ k+1
where Uik+1 and Ui+1 denote the intermediate numerical solution of u at time k+1 on spatial
¯
¯
point i and i + 1, respectively, while µ(U k+1 ) is the total variation measure of U k+1 . As
¯
¯
a result, by checking whether the increment of high frequency ∆µ(U k ) = µ(U k+1 )
¯
µ(U k )

−

exceeds η or not, we can decide whether to apply the PDE transform.

Equation (2.20) in the second step is given by the PDE transform to suppress Gibbs’
oscillations. We implement the PDE transform in frequency domain, which results in

ˆ
ˆˆ
¯
U k+1 = LU k+1 ,

(2.22)

ˆ
ˆ
¯
¯
ˆ
where U k+1 and U k+1 are the frequency responses of U k+1 and U k+1 respectively, while L
is the kernel of the fast PDE transform in the Fourier domain. In this work, we adopt the
simple implementation of PDE transforms of dm = 1 and di = 0, (i = 1, · · · , m − 1) as well
as = 0, as mentioned before. As a result, the design of the PDE transform depends on the
highest order of l = 2m and propagation time τ , whose properties have been analyzed in the
last section.
It should be noted that the FPM is normally associated with periodic boundary conditions. However, most conservation law equations admit non-periodic boundary conditions. We
handle non-periodic computation domains by symmetrical extension to an auxiliary domain.

50

This approach works well as shown in our previous work [151]. Speciﬁcally, the symmetrical
extension of the computation domain in 1D is discussed here. Suppose the original domain
is discretized by N + 1 grid points indicated by i = 1, · · · , N + 1, then the computation
domain is doubled with 2N grid points and f (u) and u are extended by

{f (u)}i = {f (u)}2N +2−i ,
{u}i = {u}2N +2−i ,

i = N + 2, · · · , 2N,
i = N + 2, · · · , 2N,

(2.23)
(2.24)

where {·}i corresponds to the value at the ith pixel. Extensions in 2D and 3D are straightforward as shown Eqs. (2.49)-(2.54) in Section 2.3.3.
After the entire time integration has been completed, we adopt a cosmetic post-processing
ﬁlter as introduced by Gottlieb et al. [76] to make the solution more presentable. In our
numerical experiments, we employ the Lagrange-4 as the post processing ﬁlter. The reader
is referred to our earlier work [152] for its implementation detail.

2.3

Numerical test and validation

The performance of the proposed PDE transform on the solution of hyperbolic conservation
law systems is validated through test examples in this section. A number of standard linear
and nonlinear benchmark problems are studied in the present work, from scalar conservation law systems including the linear advection equation, Burgers’ equation and problem
with non-convex ﬂux, to one dimensional (1D) and two dimensional (2D) Euler equations,
including shock tube (Sod’s and Lax’s) problems, 1D and 2D shock-entropy interaction, as
well as 2D shock-vortex interaction.

51

The proposed PDE transform can be applied directly on those problem with periodic boundary conditions. However, for non-periodic boundary conditions, the computation
domain is symmetrically doubled to convert to a periodic one, as in Example 3 (Inviscid
Burgers’ equation), Example 4 (Non-convex ﬂux), Example 5 (Sod’s and Lax’s problems),
Example 6 (1D Shock-entropy interaction), Example 7 (Shu-Osher’s problem), Example 8
(2D Shock-entropy interaction) and Example 9 (2D Shock-vortex interaction).
In our implementation, the wavenumber w is set by wq = ( 2πq )/∆, q = − N , · · · , 1, · · · , N −
2
2
N
1 where L is the length of computation domain, N is the number of discretized grid points
and ∆ = L/N is the grid spacing. This series of wavenumbers ranges over [−π/∆, π/∆). The
design of PDE transform depends on the order n and propagation time τ . Table 2.1 lists the
order l and propagation time τ used in each test problem. To better interpret the frequency
response, we rescaled the wavenumber w from [−π/∆, π/∆) to [−π, π). Consequently, the
propagation time τ is rescaled to τ ∗ = τ /(∆)l correspondingly.
It is noted that the appropriate selection of order l = 2m and propagation time τ of the
PDE transform is quite subtle. On the one hand, the order should not be too high, while the
propagation time should not be too small. Otherwise, the PDE transform has little eﬀect in
suppressing the Gibbs’ oscillation of the solution. On the other hand, the order should not
be too low, while the propagation time should not be too large. Otherwise, they result in
too much dissipation. Owing to diversity of hyperbolic conservation law systems, the best
choices of order and propagation time are problem dependent, which are listed in the table
2.1.
We would like to point out that both the PDE transform and hyperbolic conservation law
systems are time dependent. To avoid possible confusion, we denote “τ ” as the propagation
time of the PDE transform and designate “t” as evolution time of the conservation law
52

Table 2.1: The order and propagation time of the PDE transform used in test examples.
Example No

1

Case

Order

Propagation time
π π
when w ∈ [− ∆ , ∆ )

κ = 5, 10
κ = 20, 25

6
12
2
4
6
8
6
6
4
6
6
6
12
12
6
40
12
10
2
10
6
10

1.0E-13
1.0E-20
1.0E-05
1.0E-12
1.0E-16
1.0E-20
3.0E-15
7.5E-10
5.5E-14
6.0E-18
1.0E-09
1.0E-10
3.0E-27
1.5E-30
5.0E-16
1.0E-95
2.0E-34
1.0E-23
5.0E-06
3.5E-23
1.0E-14
1.0E-24

2

3

ul = 1, ur = 0
ul = 0, ur = 1

4
5
6

Sod’s
Lax’s
κ = 18
κ = 32

κ = 60
7
8
9

Rescaled
propagation time
when w ∈ [−π, π)
6.8E-03
4.7E-03
1.6E-01
2.6E-04
4.4E-04
7.2E-04
1.2E-03
5.3E-02
1.0E-03
5.9E-05
4.4E-03
4.4E-04
3.4E-06
7.0E-06
1.0E-03
1.7E-13
3.8E-06
1.1E-05
2.0E-02
1.2E-05
3.3E-04
1.1E-03

system in the rest of the paper.

2.3.1

Scalar conservation law systems

We ﬁrst consider 1D scalar conservation law systems, whose governing equation is expressed
in the form of

ut + f (u)x = 0,

53

(2.25)

1
with f (u) as a function of u. Three types of f (u), including f (u) = u and convex ﬂux 2 u2
1
as well as non-convex ﬂux 4 (u2 − 1)(u2 − 4), are studied in the present work.

2.3.1.1

Example 1 (Linear advection equation with Sine-Gaussian wavepacket)

. The ﬁrst test example is the linear advection equation given by

ut + cux = 0,

−1 < x < 1,

u(x, 0)

periodic,

(2.26)

= u0 (x),

where u0 (x) is the initial value. We set u0 (x) as [206]

u0 (x) = sin[2πκ(x − x0 )]e

−

(x−x0 )2
2σ 2

(2.27)

where the parameter of wavenumber κ can be tuned to produce highly oscillatory wavepacket,
x0 is the initial location of wavepacket center and σ is a constant regularizing the width of the
√
wave packet. We set x0 = 0 and σ = 2/10 such that the tail of wavepacket is constrained
over the computation domain [−1, 1]. It is straightforward to show that the exact solution
is given by u(x, t) = u0 (x − t), which is a translation of the initial solution at a unit speed,
resulting in the fact that the initial wavepacket repeats itself every two time units over the
computation domain [−1, 1]. The mesh size is chosen as N = 128 to implement the PDE
transform while the time increment is selected as small as 10−4 to ensure that the time
discretization error is negligible.
We vary wavenumber κ to explore the accuracy and stability of the PDE transform strategy for the wavepacket propagation. For low frequency waves, a low-order method can work
well to obtain accurate solutions. However, as the wavenumber κ increases, the wavepacket
54

1

Numerical
Exact

0.5

0

−0.5

−1
−1

−0.5

0

0.5

1

(a) κ = 5, N = 128
1

Numerical
Exact

0.5

0

−0.5

−1
−1

−0.5

0

0.5

1

(b) κ = 10, N = 128
Figure 2.1: Results for κ = 5, 10 from the PDE transform for the advection equation with
the Sine-Gaussian wave packet (t=100, ∆t = 10−4 ).

55

1

Numerical
Exact

0.5

0

−0.5

−1
−1

−0.5

0

0.5

1

(c) κ = 20, N = 128, before the interpolation
1

Numerical
Exact

0.5

0

−0.5

−1
−1

−0.5

0

0.5

1

(d) κ = 20, N = 128, after the interpolation
Figure 2.2: Results for κ = 20 from the PDE transform for the advection equation with the
Sine-Gaussian wave packet (t=100, ∆t = 10−4 ).

56

Table 2.2: L∞ error for long-time integration of the Sine-Gaussian wavepacket.
Time
PDE transform

κ = 20
κ = 25
CFOR-Hermite κ = 20
κ = 25

10
2.69E-7
8.10E-7
2.78E-7
1.51E-6

20
5.38E-7
1.62E-6
5.56E-7
3.02E-6

50
1.34E-6
4.06E-6
1.39E-6
7.55E-6

80
2.15E-6
6.50E-6
2.22E-6
1.21E-5

100
2.67E-6
8.47E-6
2.78E-6
1.51E-5

becomes more and more oscillatory. It is challenging for low resolution schemes to resolve
the highly oscillatory advective wavepacket. Therefore, one needs to resort to a high order
method. In our test, we apply the 6th-order PDE transform for low frequency wavepackets, i.e., κ = 5 and 10, and 12th-order PDE transform for high frequency wavepackets, i.e.,
κ = 20 and 25.

Figures 2.1(a), (b) and 2.2(a) show the computational results and exact solutions at
time t = 100 for three diﬀerent frequencies κ = 5, 10 and 20. It is obvious that there is
no visual diﬀerence between the numerical and exact solutions in Figures 2.1(a) and (b).
However, for the case of of κ = 20, the solution does not show the peaks and valleys of
the exact wavepacket. One may regard the solution as being imperfect. In fact, this is not
true. The L∞ error listed in Table 2.2, is as small as 2.67E-6. Therefore the numerical
solution is extremely accurate. The truth is that, because so few grid points are used
in our computation, there are not enough grid points to fully represent the wavepacket.
Nevertheless, the origin information of the wavepacket is perfectly built in the solution. To
illustrate this point, we interpolate the present solution to a denser grid (N = 256) as shown
in Figure 2.2(b) by using our DSC algorithm [167, 177, 183]. Indeed, we are able to fully
restore all extrema in the wavepacket based on the information presented in Figure 2.2(a).
This conﬁrms that our scheme is still able to perform well for high frequency waves under a

57

very small ratio of grid points over wavelength.
It is also meaningful to check the long time integration of the numerical method. Table 2.2
lists the numerical L∞ -error from t = 10 to t = 100. One can tell that the numerical errors
are under good control during the time integration. As a comparison, we also list in Table
2.2 the results of a CFOR-Hermite method proposed in our earlier work [206]. The CFORHermite method is based on the local spectral wavelet approach for the spatial discretization
and the conjugate ﬁlter method for oscillation reduction. It has been intensively validated
for solving the Navier-Stokes equation and integrating conservation law systems [206].
In the work of [206], the CFOR-Hermite method employs 100 points but takes the same
time step ∆t = 10−4 . From Table 2.2, one can notice that the errors from the PDE transform
is slightly smaller than those from the CFOR-Hermite method. We are quite conﬁdent that
the results listed in Table 2.2 are some of the best when κ = 25. It can be concluded that the
numerical accuracy and long time stability are well resolved by the proposed PDE transform.

2.3.1.2

Example 2 (Linear advection equation with wave combination)

Another example is also a linear advection equation given by [95]

ut + ux = 0,

−1 < x < 1,

u(x, 0)

periodic,

(2.28)

= u0 (x),

58

where u0 (x) is the initial value given by

u0 (x) =


 1
 (G(x, β, z − δ) + G(x, β, z + δ) + 4G(x, β, z)),
 6






 1,




 1 − |10(x − 0.1)|,




 1 (H(x, α, a − δ) + H(x, α, a + δ) + 4H(x, α, a)),
 6






 0,

−0.8 ≤ x ≤ −0.6,
−0.4 ≤ x ≤ −0.2,
0 ≤ x ≤ 0.2,

(2.29)

0.4 ≤ x ≤ 0.6,
otherwise.

The functions G and H are set as
2

G(x, β, z) = e−β(x−z) ,
(2.30)
H(x, α, a) =

max(1 − α2 (x − a)2 , 0),

where z = −0.7, δ = 0.005, β = log2 , a = 0.5, and α = 10.
2
36δ

The initial value of this problem is a smooth but narrow combination of a Gaussian,
a square wave, a sharp triangle wave and a half ellipse. It is easy to show that the exact
solution is given by u(x, t) = u0 (x − t), which is a translation of the initial solution at a
unit speed. It is well known that due to the contact discontinuity, the propagation by the
linear advection equation leads to unphysical Gibbs’ oscillations which may be induced by
exponentially small numerical errors and their subsequent ampliﬁcation. In Figure 2.3, the
numerical results obtained by using the PDE transform of second, fourth, sixth and eighth
orders are demonstrated at t = 8 with 256 grid points. It is interesting to observe that for
this simple problem, not only the second order PDE transform works well, but also many
higher order PDE transforms produce satisfactory results.

59

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0
−1

−0.5

0

0.5

0
−1

1

(a) The 2nd-order

−0.5

0

0.5

1

(b) The 4th-order

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0
−1

−0.5

0

0.5

0
−1

1

(c) The 6th-order

−0.5

0

0.5

1

(d) The 8th-order

Figure 2.3: Results from the PDE transform of various orders for the advection equation (
t = 8, ∆t = 0.001, 256 grid points).

60

1
0.8
0.6
0.4
0.2
Numerical
Exact

0
−1

−0.5

0

0.5

1

1.5

2

(a)

1

Numerical
Exact

0.8
0.6
0.4
0.2
0
−1

0

1

2

3

(b)
Figure 2.4: Results from the 6th-order PDE transform for the inviscid Burgers’ equation at
time t = 2. (a) ul = 1, ur = 0, 257 grid points; (b) ul = 0, ur = 1, 129 grid points;

61

2.3.1.3

Example 3 (Inviscid Burgers’ equation)

In this example, we test the performance of the present PDE transform for the most classical
1
model, inviscid Burgers’ equation with convex ﬂux f (u) = 2 u2 and Riemann type of initial

conditions with constant ul on the left and ur on the right.
(3a) First we consider the Riemann initial value with ul > ur as


 ul = 1,
u(x, 0) =

 u = 0,
r

x ≤ 0,
(2.31)
x > 0.

This problem has been studied by numerous researchers because it is a standard benchmark
problem in hyperbolic conservation laws. The exact solution is given by a shock wave with
a constant velocity σ, i.e.

x
t

< σ,


 u ,
r

u(x, t) =



 ul ,

x
t

> σ,

(2.32)

where
σ=

f (ul ) − f (ur )
1
= .
ul − ur
2

It is observed that this problem has a non-periodic boundary condition. Consequently,
the computational domain needs to be symmetrically doubled in the x-direction to obtain
periodic boundary condition when applying the Fourier pseudospectral method [151] pairing
with the PDE transform. Since the initial condition is piecewise constant and exact solution
does not involve large oscillation, a low-order scheme is suitable to render a satisfactory
numerical solution. In this case, we apply second, fourth and sixth order PDE transforms
to resolve the problem. It is found that all the three schemes work well. The numerical
62

results from the sixth order PDE transform are plotted in Figure 2.4(a) at time t = 2. It is
demonstrated that the Gibbs’ oscillation is well resolved and the shock front, which moves to
x = 1, is well captured. Although this problem prefers low order shock-capturing methods
[95], it is concluded that the PDE transform based FPM method works well for this problem
too.
(3b) As another example, we check the present method by using the Riemann type initial
value with ul < ur as




ul = 0,

x < 0,


 u = 1,
r

u(x, 0) =

x ≥ 0.

(2.33)

The exact solution is given by a rarefaction wave



 0,



x
u(x, t) =
 G( t ),




 1,

x
t

< f (ul ),

f (ul ) < x < f (ur ),
t

(2.34)

f (ur ) < x .
t

with f (ul ) = ul , f (ur ) = ur and G( x ) = f
t

−1 x
(t)

= x.
t

Similar to Case (3a), the computational domain is symmetrically doubled in the xdirection to obtain periodic boundary condition before applying the PDE transform and
FPM. Furthermore, the piecewise constant initial condition and exact solution imply that
this problem can be well resolved by a low-order shock-capturing scheme. In our approach,
we apply the sixth order PDE transform, paired with the FPM, to solve the equation. The
numerical results from the sixth order PDE transform are plotted in Figure 2.4(b) at time
t = 2. From the ﬁgure, it is seen that the rarefaction fan over [0,2] is free of oscillation. The

63

feature of the solution is well resolved.

2.3.1.4

Example 4 (Non-convex ﬂux)

We next consider a problem with a non-convex ﬂux to test the convergence to the physically
correct solution. The non-convex ﬂux is given by

1
f (u) = (u2 − 1)(u2 − 4)
4

(2.35)

with a Riemann initial condition


 ul = −3,
u(x, 0) =

 u = 3,
r

x < 0,
(2.36)
x ≥ 0.

The exact solution is given by

u(x, t) =



 u,
 l






 G( x ),
t

x
t

≤ f (ul ),

f (ul ) < x < 0,
t


 −G(− x ),


t





 ur ,

0<

x
t

(2.37)

< f (ur ),

f (ur ) ≤ x ,
t

with f (ul ) = −19.5, f (ur ) = 19.5 and G(u) is the solution of f (G(u)) = u in the convex
part of f , which is | u |>
√
x = 0 and G(0) = 2.5.

5/6. It should be noted that the solution is discontinuous at

The exact solution and more detailed information about this problem are given in [90].
This problem is relatively more complicated than the convex ﬂux case. In the literature,
commonly reported numerical result is at t = 0.04 [95]. The numerical results of the PDE

64

3
2

Numerical
Exact

1
0
−1
−2
−3
−1

−0.5

0

0.5

1

(a) The 4th order PDE transform (129 grid points)

3
2

Numerical
Exact

1
0
−1
−2
−3
−1

−0.5

0

0.5

1

(b) FPM-RSK (129 grid points)
Figure 2.5: Comparison between numerical results from the PDE transform and FPM-RSK
method for the problem with a non-convex ﬂux ( 129 grid points, t = 0.04, ∆t = 0.0005).

65

3
2

Numerical
Exact

1
0
−1
−2
−3
−1

−0.5

0

0.5

1

(a) The 6th order PDE transform (65 grid points)
3
2

Numerical
Exact

1
0
−1
−2
−3
−1

−0.5

0

0.5

1

(b) FPM-RSK (65 grid points)
Figure 2.6: Comparison between numerical results from the PDE transform and FPM-RSK
method for the problem with a non-convex ﬂux ( 65 grid points, t = 0.04, ∆t = 0.0005).

66

transform are showed in Figures 2.5(a) and 2.6(a). It is seen that the shock front at x = 0
is almost exactly captured. Particularly, there is no numerical solution point located on the
discontinuity.
As a comparison, we consider the Fourier pseudospectral method (FPM) with the regularized Shannon kernel (RSK) approach, which was extensively validated in our earlier work
[151]. In Figures 2.5(b) and 2.6(b), this approach is labeled as FPM-RSK. It is obvious
that compared with the FPM-RSK, the present PDE transform yields a more satisfactory
resolution. We believe that our results from the PDE transform based FPM approach are
some of the best ever reported for this classic problem.

2.3.2

1D Euler systems

In this subsection, we carry out numerical experiments by using the proposed PDE transform
scheme for the 1D Euler equation of gas dynamics. In one dimension, the Euler equation
takes the form [95, 151]

Ut + F (U )x = 0

(2.38)

with




 ρ 




U =  ρu  ,




E









 ρ 
 0 








F (U ) =  ρu  u +  p  ,








E
pu

67

(2.39)

where ρ, u, p and E denote the density, velocity, pressure and total energy per unit mass

1
E = ρ[e + u2 ],
2

(2.40)

respectively. Here, e is the speciﬁc internal energy. For an ideal gas with the constant speciﬁc
heat ratio (γ = 1.4) considered here, one has

e =

p
.
(γ − 1)ρ

(2.41)

In the following, we consider two well-known Riemann problems.

2.3.2.1

Example 5 (Sod’s and Lax’s problems)

Here we apply the PDE transform based FPM on two shock tube problems, i.e., Sod’s
problem and Lax’s problem, which are both standard benchmark tests. In fact, due to their
simple proﬁles, these problems favor low order schemes.
Sod’s problem is a special case of shock tube problem with velocities on both sides of
the discontinuity being set to zero. It is often used as a test case for validation of numerical
shock capturing schemes, because analytical solutions are available. The initial and boundary
conditions in the computation domain and boundary for Sod’s problem are given by

(ρ, u, p)t=0 =



 (1,

0,


 (0.125,

0,

1),

−5 ≤ x < 0,
(2.42)

68

0.1),

0 ≤ x ≤ 5;

1.2

1.2

Numerical
Exact

1

1
Pressure

0.8
Density

Numerical
Exact

0.6

0.8
0.6

0.4

0.4

0.2

0.2

0
−5

0
Position

0
−5

5

(a)

0
Position

5

(b)

Figure 2.7: Numerical results from the 6th-order PDE transform for Sod’s problem (t = 1.5,
∆t = 0.02, 129 grid points). (a) Density; (b) Pressure.
while those for Lax’s problem are given by

(ρ, u, p)t=0 =



 (0.445,

 (0.5,

0.698,

3.528),

−5 ≤ x < 0,
(2.43)

0,

0.571),

0 ≤ x ≤ 5.

Although there are coupled equations in the system, which make it more complicated
than scalar equations, the Sod and Lax problems involve multiple piecewise constant solution
without much oscillation. Therefore, these problems favor low-order numerical methods. In
our study, we utilize the 6th-order PDE transform method with the FPM to integrate these
shock tube equations. Our results of density and pressure for Sod’s problem are depicted
in Figure 2.7. It is perceived that three characteristics of Sod’s problem, including the
rarefaction wave, the contact discontinuity and the shock discontinuity are well resolved and
captured.
Figures 2.8(a) and (b) illustrate the present numerical results for Lax’s problem. Similar
69

1.4
1.2

4

Numerical
Exact

3
Pressure

1
Density

Numerical
Exact

0.8
0.6
0.4

2

1

0.2
0
−5

0
Position

0
−5

5

(a)
1.4
1.2

5

(b)
4

Numerical
Exact

Numerical
Exact

3
Pressure

1
Density

0
Position

0.8
0.6
0.4

2

1

0.2
0
−5

0
Position

0
−5

5

(c)

0
Position

5

(d)

Figure 2.8: Comparison of numerical results from the 6th-order PDE transform and the
FPM-RSK for Lax’s problem (t = 1.5, ∆t = 0.02, 129 grid points). (a) Density from the
PDE transform; (b) Pressure from the PDE transform; (c) Density from the FPM-RSK; (d)
Pressure from the FPM-RSK.

70

to the Sod’s case, three characteristics of Lax’s problem, including the rarefaction wave, the
contact discontinuity and the shock discontinuity are also well captured. It can be seen that
the present method gives a good resolution in both cases.
As a comparison, we plot the results obtained by the FPM-RSK approach [151] in Figures 2.8(c) and (d). As discussed earlier, the FPM-RSK scheme is a robust approach for
hyperbolic conservation laws. It is seen from Figure 2.8 that the present PDE transform
performs at least as good as, if it is not better than, the FPM-RSK approach for shock tube
problems.
All of the preceding test examples, except for the high-frequency case in Example 1, favor
low order shock capturing schemes. In such a case, although the present method, as well as
other high-order shock capturing schemes, works extremely well, it does not have a cutting
edge advantage over low order methods. In the next three examples, we consider a class
of problems that require high-order shock capturing methods methods for eﬃcient spatial
discretizations. Whereas, it will be extremely diﬃcult, if it is not impossible, for low order
methods to resolve this class of problems.

2.3.2.2

Example 6 (1D Shock-entropy interaction)

The interaction between a Mach 3 right-moving shock and an entropy wave of small amplitude in a one-dimensional ﬂow, which is a standard test problem [95, 151, 85], is investigated
in this example. This problem is important because of its relevance to the interaction of shock
and turbulence. We take the computation domain over [0,9] and the initial and boundary

71

0.1

0.05

0

−0.05

−0.1
4

5

6

7

8

9

7

8

9

7

8

9

(a) κ = 18
0.1

0.05

0

−0.05

−0.1
4

5

6

(b) κ = 32
0.1

0.05

0

−0.05

−0.1
4

5

6

(c) κ = 60
Figure 2.9: The amplitude of entropy waves in the 1D shock-entropy interaction problem for
diﬀerent frequencies resulting from the 12th-order PDE transform.

72

0.1

0.05

0

−0.05

−0.1
4

5

6

7

8

9

7

8

9

(a)
0.1

0.05

0

−0.05

−0.1
4

5

6

(b)
Figure 2.10: The amplitude of entropy waves in the 1D shock-entropy interaction problem for
κ = 32 resulting from PDE transforms of low and extremely high-order. (a) The 6th-order
PDE transform; (b) The 40th-order PDE transform.

73

conditions in the computation domain and boundary are given by

2.629369,

 − sin(κx)
 (e
,

(ρ, u, p)t=0 =



 (3.85714,

0,

10.33333),

0 ≤ x < 0.5,
(2.44)

1.0),

0.5 ≤ x ≤ 9.

The parameters and κ are the amplitude and the wave number of the entropy wave before
the shock. In this numerical experiment, the small amplitude of the pre-shock entropy is
kept the same, i.e.

= 0.01 while the wave number κ is varied. The amplitude of the

ampliﬁed wave after the shock can be given by linear analysis [123], which is a constant of
0.08690716.
This problem of shock and entropy interaction becomes more and more challenging as the
frequency κ increases. The diﬃculty lies in the fact that it is hard to distinguish the ampliﬁed
high frequency entropy wave from the spurious oscillations. Low order numerical schemes
may dramatically damp the amplitude of the transmitted high frequency wave. Even some
high order schemes encounter the same numerical diﬃculty. A satisfactory numerical scheme
should possess the quality of suppressing Gibbs’ oscillation while preserving amplitude of the
entropy wave as well as capturing the shock [151, 95, 85].
In our numerical test, the shock is set to move from x = 0.5 to x = 8.5. For the
convenience of a comparison with previous results, we only display the results at interval
[3.0, 9.0]. Furthermore, for the purpose of discharge transient waves due to the non-numerical
initial shock proﬁle, we plot the length of the ampliﬁed entropy waves in the same manner
as that in Ref. [95]. Additionally, we would like to point out two nontrivial details in
the illustration. The ﬁrst one is that the plotted results are obtained by an interpolation
of our ﬁnal numerical results from a coarse grid to a denser grid, which has already been
adopted for plotting the high-frequency Sine-Gaussian wavepacket in Example 1. This is
74

necessary because our computational grid is too coarse to show all the peaks and valleys of
our numerical results. The second one is that a post processing ﬁlter, which was discussed
in Section 2.2.3, is employed to eliminate the oscillations near the shock when we present
the ﬁnal result after the completion of the entire time integration.
First, we test the shock capture for a case with relatively low frequency κ = 18. We
set 513 grid points over the computation domain [0, 9] to implement the PDE transform
coupled with the Fourier pseudospectral method. It is diﬃcult for a low-resolution method,
such as a ﬁrst order or a second order shock capturing method to preserve the amplitude of
high-frequency entropy wave and suppress oscillation [95].
In our numerical test study, we found that the low-order PDE transform does not work
well for this problems, although the FPM is employed. Low order PDE transforms either
damp the amplitude of the high-frequency entropy wave or cannot suppress Gibbs’ oscillations. In contrast, we found that a moderately high order PDE transforms with suitable propagation time perform much better. Figure 2.9(a) shows the amplitude of the shock entropy
of κ = 18 obtained by using the 12th-order PDE transform. It can be seen that the entropy
waves almost fully span two strips bounded by two analytical solid lines y = 0.08690716 and
y = −0.08690716, which indicates that the amplitude of ampliﬁed wave is well preserved. As
the ampliﬁed entropy wave is monochromic post the shock, it is appropriate to characterize
the resolution of the present method by points per wavelength (PPW). A further simple
calculation on the ampliﬁed entropy wave tells that the resolution is about 5.3 PPW, which
is among the best results for this shock entropy problem, to our best knowledge.
Next we increase the frequency to κ = 32. As there are many more ampliﬁed entropy
wave post the shock, it is diﬃcult to use 513 grid points to maintain the amplitude of highfrequency wave post shock. Instead, we deploy a mesh of 1025 grid points so as to render a
75

better result. In this case, we still employ a 12th-order PDE transform but with a diﬀerent
propagation time from κ = 18. One still needs to be careful to design an appropriate PDE
transform for this more complicated problem. Figure 2.9(b) shows the amplitude of the shock
entropy wave of κ = 32 by a 12th-order PDE transform. It is perceived that the entropy
wave is still able to almost fully span two strips bounded by two solid lines y = 0.08690716
and y = −0.08690716. Furthermore, the waves at shock front are free of oscillation.
Finally, we consider a more challenging situation by setting a larger κ. We use 2049 grid
points in our spatial discretization of the computation domain [0, 9]. Figure 2.9(c) presents
the result of the case for κ = 60. It is clear that the waves post shock are well preserved while
the waves at shock front are free of oscillation and well maintained. It also worthwhile to
mention that the resolution for this case is about 5.9 PPW, which is still the best record at
such a high κ value to our knowledge. Our results indicate that the present PDE transform
is able to distinguish the high-frequency entropy waves from spurious oscillation. The above
test cases of diﬀerent frequencies imply that the PDE transform can be a powerful tool
for solving hyperbolic conservation law systems involving the interaction of shock and high
frequency entropy waves.
As stated earlier, low order shock capturing schemes will encounter diﬃculties for the
above test cases. In fact, not all PDE transforms are suitable for these problems. The
inappropriate selections of the PDE transform will either pollute the amplitude of entropy
waves or cannot suppress the oscillation. We illustrate this point by using the setting of
κ = 32 and 1025 grid points. Figure 2.10(a) demonstrates the result from a 6th-order PDE
transform, which illustrates the eﬀect of a low order PDE transform. Although the waves in
the shock front is perfectly kept, the amplitude of the entropy waves post shock is severely
damped. We found that second order and fourth order PDE transforms perform even worse
76

4

3

3

Density

5

4
Density

5

2
1
0
−1

1

Numerical
Exact
−0.5

2

0
Position

0.5

0
−1

1

(a)

Numerical
Exact
−0.5

0
Position

0.5

1

(b)

Figure 2.11: Solution to the Shu-Osher problem (t = 0.47, ∆t = 0.001, 129 grid points).
(a) Obtained with the 2nd-order PDE transform; (b) Obtained with the 10th-order PDE
transform.
for this problem. As another check, we would like to investigate the performance of a very
high order PDE transform for this problem. Although a high order PDE transform works
well in preserving the high-frequency entropy waves post the shock, it may also pollute the
relatively low frequency waves in the shock front. Figure 2.10(b) shows the result from a
40th-order PDE transform. It is clear that the entropy waves in the shock front are polluted
by Gibbs’ oscillations.

2.3.2.3

Example 7 (Shu-Osher’s problem)

We now examine the performance of the PDE transform on the problem of Shu and Osher,
which is another typical case to test the capability of a numerical method in predicting
shock/entropy interactions. This problem is also initialized by a Mach 3 right-moving shock
and an entropy wave, which is in the sine form. Compared with the last problem, Shu and
Osher problem involves many diﬀerent frequencies. We consider the computational domain

77

of [−1, 1] and the initial and boundary conditions in the computation domain and boundary
are set as

(ρ, u, p)t=0 =





(3.85714,


 (1.0 + sin(κπx),
where

2.629369,

10.33333),

−1 ≤ x < −0.8,
(2.45)

0,

1.0),

−0.8 ≤ x ≤ 1,

= 0.2 is the amplitude and κ = 5 the angular wave number. Although there is no

exact solution to this problem, we took the accurate solution as the result given by the highly
accurate FPM-RSK method [151], which was calibrated by the 5th order WENO scheme [95]
for this problem.
Because of the complicated ﬂuctuation post the shock, it is diﬃcult for low order schemes
to capture various frequencies of wave post shock as well as the sharp shock front [95]. High
order methods, such as the 5th order WENO scheme [95] and the FPM-RSK method perform
extremely well [151].
In the present work, we are interested in the understanding of the performance of the
PDE transform of diﬀerent orders. In fact, through our numerical tests, we found that
when one uses the second order PDE transform, one cannot obtain a satisfactory numerical
result, no matter what propagation time is chosen. With various choices of the propagation
time, the second order PDE transform either damps the amplitude of high frequency waves
too much or cannot prevent the solution from blowing up as time evolves. Figure 2.11(a)
shows the result of a 2nd-order PDE transform. It is seen that it not only smooths the low
frequency waves but also damps the high frequency components.
By contrast, high order PDE transforms work much better. As an illustration, Figure
2.11(b) presents the result from a 10th order PDE transform. Obviously, the 10th-order
PDE transform is able to well capture the complicated high-frequency waves post the shock.
78

Furthermore, the shock front remains sharp in high-order PDE transforms. Although there
is some abnormal oscillation at the on-set of low frequency waves, the 10th order PDE
transform preserves the majority of the low frequency waves. This diﬃcult case further
validates the power of the present PDE transform.

2.3.3

2D Euler systems

In two spatial dimensions, we consider the Euler equation in the conservative form of

Ut + F (U )x + G(U )y = 0

(2.46)

with




 ρ 




 ρu 


U =
,


 ρv 




E






 ρ 







 ρu 



F (U ) = 
u + 




 ρv 






E







0 


p 

,

0 


pu






 ρ 







 ρu 



G(U ) = 
v + 




 ρv 






E


0 


0 

 , (2.47)

p 


pv

where (u, v) is the ﬂuid velocity and p is given by

1
p = (γ − 1)[E − ρ(u2 + v 2 )].
2

(2.48)

We discretize the conservative quantities (ρ, ρu, ρv, p) on the mesh and applied the PDE
transform during the time integration.
Like the one dimensional problems, the non-periodic domain in 2D can also be handled
by symmetrical extension to the auxiliary domain. Suppose the original domain is discretized

79

by N1 · N2 grid points indicated by (i, j) with i = 1, · · · , N1 + 1 and j = 1, · · · , N2 + 1. If
the 2D domain is non-periodic in only one of the two directions, suppose in x-direction, then
the computation domain is doubled with 2N1 · (N2 + 1) grid points and F (U ), G(U ) and u
are extended by

{F (U )}i,j,k = {F (U )}2N1 +2−i,j,k ,

i = N1 + 2, · · · , 2N1 ,

(2.49)

{G(U )}i,j,k = {G(U )}2N1 +2−i,j,k ,

i = N1 + 2, · · · , 2N1 ,

(2.50)

{U }i,j,k = {U }2N1 +2−i,j,k ,

i = N1 + 2, · · · , 2N1 .

(2.51)

If the 2D domain is non-periodic in both two directions, then the computation domain is
doubled with 2N1 · 2N2 grid points and F (U ), G(U ) and U are extended by

{F (U )}i,j,k =

{G(U )}i,j,k =

{U }i,j,k =



 {F (U )}2N +2−i,j,k ,
1

 {F (U )}
i,2N2 +2−j,k ,


 {G(U )}2N +2−i,j,k ,
1

(2.52)

j = N2 + 2, · · · , 2N2 ,
i = N1 + 2, · · · , 2N1 ,


 {G(U )}
j = N2 + 2, · · · , 2N2 ,
i,2N2 +2−j,k ,


 {U }2N +2−i,j,k ,
i = N1 + 2, · · · , 2N1 ,
1

 {U }
i,2N2 +2−j,k ,

2.3.3.1

i = N1 + 2, · · · , 2N1 ,

(2.53)

(2.54)

j = N2 + 2, · · · , 2N2 .

Example 8 (2D Shock-entropy interaction)

Having tested the performance of the PDE transform for 1D shock-entropy interactions, we
now consider the case in a 2D setting. The weak entropy wave makes an angle θ ∈ (0, π/2)
against the x-axis. If θ = 0, then this 2D problem essentially degenerates into the 1D shock80

0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
7.4

Linear analysis
The 10th−order PDE transform
The 6th−order PDE transform
7.6

7.8

8

8.2

8.4

Figure 2.12: The amplitude of entropy waves in the 2D shock-entropy interaction problem
obtained from PDE transforms of diﬀerent orders.

81

entropy problem studied earlier. Since now the entropy waves are oblique to the shock, this
2D problem is more challenging to resolve. Our objective is to further examine the capability
of the PDE transform for higher dimensional problems.
Given the right state of the shock as (ρr , ur , vr , pr ) = (1, 0, 0, 1), a weak entropy is added
by changing the density on the right. Here ρr is modiﬁed by

ρ = ρr e− sin(κz(θ))/pr ,

where z(θ) is related to the angle θ by z(θ) = x cos θ + y sin θ and

(2.55)

and κ are the amplitude

and wave number of the entropy wave, respectively. In this test, we choose parameters
θ = π/6,

= 0.1 and κ = 15.

For the computation domain, on the one hand, in order to implement the periodic bound2π
ary condition along the y-direction, the computation domain in y is set to be [0, κ sin θ ] pro-

vided θ = 0. On the other hand, the computation domain in x is set as [0, 9].

Since the

boundary condition in the x-direction is non-periodic, we extend the computation domain
2π
in the x-direction symmetrically from [0, 9] to [0, 18] . We deploy 32 points in [0, κ sin θ ] and

1024 points over [0, 18], which implies 513 points over [0, 9]. In our test, the shock starts at
x = 0.5 and moves up to x = 8.5. As we mentioned in the 1D case, a good numerical scheme
should be able to preserve the amplitude of the high-frequency entropy waves. It needs to
be mentioned that low-order methods do not work well for this problem [95]. Low-order
methods would severely damp the amplitude as pointed out in Ref. [95].
For the PDE transform approach, methods of diﬀerent PDE orders exhibit sharply diﬀerent behaviors. Figure 2.12 shows the performance of the 6th-order PDE transform by
checking the maximum amplitude of the ampliﬁed waves along the y-direction for grid points

82

x ∈ [7.4, 8.4], and furthermore comparing with the amplitude predicted by the linear analysis, i.e., 0.08744786. It is seen that the 6th-order method has an obvious amplitude loss.
In fact, to preserve the amplitude of the entropy waves, a high order PDE transform is
preferred. We improve the performance by resorting to the 10th-order PDE transform. The
result from the 10th-order PDE transform is also illustrated in Figure 2.12. From this result, it is obvious that the amplitude of entropy waves is well maintained in the post-shock
region. Although there are small peaks and valleys, such a trivial deﬁciency appears in other
methods as well [95, 151] and is acceptable. This test further validates the capability of the
present PDE transform for resolving the shock entropy problem in a higher dimension.

2.3.3.2

Example 9 (2D Shock-vortex interaction)

Finally, we consider the problem of the interaction between a stationary shock and vortex. This problem has attracted much attention from numerous researchers because of its
potential applications. This is another problem that poses challenges to low order numerical schemes. High order methods, such as WENO methods [95], spectral methods [151]
and CFOR-Hermite scheme [206] work well for this problem. Here we further examine the
capability of proposed PDE transforms for shock-vortex interactions.
The set up of the problem is as follows. The original computational domain is set to
[0, 1] × [0, 1] with a stationary normal shock at x = 0.5 normal to x-axis. A ﬂow with
Mach number Ms = 1.1 enters at the inlet from the left and the shock is initialized as
√
(ρl , ul , vl , pl ) = (1, 1.1 γ, 0, 1) on the left state. A vortex centered at (xc , yc ) = (0.25, 0.5)
is generated by introducing a perturbation to the original velocity ﬁeld (u, v), temperature

83

1

1.25
1.2

0.8

1.15
0.6

1.1
1.05

0.4

1
0.95

0.2

0.9
0
0

0.2

0.4

0.6

0.8

1

(a) t = 0.05
1
1.25
0.8

1.2
1.15

0.6

1.1
1.05

0.4

1
0.2

0.95
0.9

0
0

0.2

0.4

0.6

0.8

1

(b) t = 0.2
Figure 2.13: The pressure proﬁle of 2D shock-vortex interaction problem from the 10th-order
PDE transform at t = 0.05, 0.2 (20 contours).

84

1

1.5

0.8

1.4

0.6

1.3

0.4

1.2

0.2

1.1

0
0

1
0.2

0.4

0.6

0.8

1

(a) t = 0.35

1

1.3
1.25

0.8

1.2
0.6

1.15
1.1

0.4

1.05
0.2

0

1
0.95
0.6

0.8

1

1.2

1.4

(b) t = 0.6
Figure 2.14: The pressure proﬁle of 2D shock-vortex interaction problem from the 10th-order
PDE transform at t = 0.35, 0.6 (20 contours).

85

T and entropy S. We denote the perturbation by (u , v , T , S ) , which is given as

(u , v ) =

2
τ eα(1−τ ) (sin θ, − cos θ)

T
S

where, τ = rrc , r =

= −

(γ − 1) 2 2α(1−τ 2 )
e
,
4γα

= 0,

(2.56)
(2.57)
(2.58)

(x − xc )2 + (y − yc )2 and θ = tan−1 [(y − yc )/(x − xc )]. By the relation

of T = p/ρ and S = ln ρp , the perturbation ρ and p to initial ρ and p can be derived from
γ
T and S . Here

and α are the strength and decay rate of the vortex and rc is a parameter

to regulate the strength of the vortex. In our test, we take = 0.3, α = 0.204 and rc = 0.05.
The computational domain [0, 1] × [0, 1] is extended to [0, 2] × [0, 2] and discretized with
an even-spacing Cartesian mesh. We use 129 even-spacing grid points in y domain of [0, 1]
and 257 grid points in x domain of [0, 1]. The mesh in the x direction is shifted by the
Robert transformation [5] to deploy more mesh points towards the stationary shock. The
upper and lower boundaries are imposed with the reﬂective boundary conditions. The time
step in our integration is regulated by the CFL condition (CFL=0.5).
This problem prefers moderately high order PDE transforms because of its complicated
feature. In fact, we cannot ﬁnd satisfactory results from low-order PDE transforms yet. As
an illustration, Figures 2.13 and 2.14 depicts our results obtained with a 10th-order PDE
transform. We plot the pressure proﬁle with 20 contours at t = 0.05, 0.2, 0.35, 0.6.
From Figures 2.13 and 2.14, it is seen that the solution is essentially free of oscillations.
It can also be observed that the deformation of the vortex at the shock and the bifurcation
after the shock are quite well captured.

86

2.4

Chapter conclusion remarks

Systems of hyperbolic conservation laws are of fundamental importance in science and engineering. The construction of accurate, eﬃcient and reliable numerical schemes has been
an active research topic in applied mathematics over the past half a century. Due to a vast
variety of complex problems with extremely diversiﬁed physical origins, many hyperbolic
conservation law systems remain a challenge to mathematical methods. Typical examples
include shock-turbulence interactions which involve a wide range of spatial scales and possible blow ups due to the successive ampliﬁcation of exponentially small numerical errors
during the time integration.
The present work introduces the use of evolution PDEs as a means to eﬃciently suppress
Gibbs’ oscillations in the numerical solution of hyperbolic conservation laws. Speciﬁcally,
during the time integration of a hyperbolic conservation law system, an intermediate numerical solution at a given time step may be used as an initial data for a special evolution PDE.
Then the solution of such an evolution PDE is accepted as an updated numerical solution
at the given time step. Our approach involves the use of the PDE transform, a technique
developed in our recent work [171, 169] for the mode decomposition of signals, image, and
data. The PDE transform is based on a family of arbitrarily high order nonlinear PDEs
originally introduced by Wei in 1999 [178] and a recursive scheme for reinitializing the input
data [170, 169]. Like the wavelet transform, the PDE transform is able to decompose signals,
images and data into functional modes, such as edges, trend, texture and feature with controllable frequency ranges and time-frequency localizations, which correspond to appropriate
multiresolution analysis in the physical domain. Similarly, the PDE transform also has a
prefect reconstruction of original signals, images and data. The PDE transform has found

87

its success in signal processing [171, 169], image analysis [171, 169] and biomolecular surface
construction [203].
It may appear computationally ineﬃcient to suppress Gibbs’ oscillations at a time step
of integrating a hyperbolic conservation law equation by solving another evolution PDE.
However, this is not true. Two techniques are proposed to improve the eﬃciency of the
present approach. First, we use an adaptive measure of total variations to automatically
determine whether the PDE transform is needed at each time step. Additionally, we utilize a
fast PDE transform, which oﬀers the analytical solution of an arbitrarily high order evolution
PDE in the Fourier representation. This technique bypasses the stability constraint of solving
high order evolution PDEs. Consequently, the present PDE transform algorithm is at least
as eﬃcient as our previous windowed Fourier pseudospectral method (FPM) [151] and is
slightly more eﬃcient than our earlier conjugate ﬁlter oscillation reduction (CFOR) scheme
[85, 181, 205, 206].
To be more speciﬁc about the eﬃciency, since the FPM with the fast Fourier transform
(FFT) is utilized as the spatial discretization, the complexity of the present PDE transform
coupled with the FFT is of O(N ln N ). This feature endows the proposed method with high
eﬃciency, which is desirable for large scale problems in scientiﬁc and engineering applications.
A variety of benchmark tests are employed in the present work to validate the proposed
approach, ranging from scalar conservation law systems to Euler equations in one and two
spatial dimensions. Among these problems, Examples 2-5 typically prefer low-order shock
capturing schemes; whereas Examples 1, and 6-9 are well known to require high order numerical methods. For example, low order schemes will severely damp the amplitude of the
entropy waves in the shock-entropy interaction described by the Euler equation. The proposed PDE transform based FPM works extremely well for these two types of problems.
88

For instance, it provides some of the best results for solving the Burgers’ equation with
non-convex ﬂux. Furthermore, only about 5 points per wavelength (PPW) is needed for the
present approach to handle the interaction of shock-entropy waves and shock-vortex interactions. To our knowledge, the only other shock-capturing schemes that have demonstrated
their ability of operating at 5 PPW are the CFOR scheme [85, 181, 205, 206] and the windowed FPM [151] proposed in our previous work. The performance of the proposed method
is compared with those of the CFOR scheme [206] and the windowed FPM [151].
In order to make the present PDE transform based FPM working well for all of the
above mentioned problems, we adjust two controlling parameters: the highest order of the
PDE transform and the propagation time. In fact, the former is the primary parameter
and the latter is less important although indispensable. It is found that the FPM behaves
as a low order shock capturing scheme when it is coupled with a relatively low order PDE
transform; while it behaves as a high order shock capturing scheme when it is coupled
with a moderately high order PDE transform. The preferred orders of the PDE transform
for these test examples are in a similar range of orders used in our PDE transform based
molecular surface construction [203]. Unlike the signal decomposition, which requires the
use of extremely high order PDE transforms [171, 169], the present systems of hyperbolic
conservation laws do not need extremely high order PDE transforms. In fact, it is found
that the use of extremely high order PDE transforms leads to unphysical oscillations. The
selection of PDE transform parameters for all test examples is summarized in the Appendix,
where the propagation time is optimized for each given order of the PDE transform.
In the present work, we employ the fast PDE transform essentially as a spectral ﬁlter
method. However, PDE transform with nonlinear geometric PDEs, such as mean curvature
ﬂows, Willmore ﬂows, surface diﬀusion ﬂows, and potential driven geometric ﬂows [14, 37, 38],
89

can be similarly used. The advantage of geometry ﬂow equation based nonlinear PDE
transform is that it enables the embedding of geometric features. It remains an open question
how geometric structures and other mathematical properties, i.e., surface decreasing and
volume preserving, will impart the solution of the systems of hyperbolic conservation laws.
The investigation of this aspect is beyond the scope of the present work, which serves only
as an introduction to the PDE transform approach for systems of hyperbolic conservation
laws.
The working principle of the PDE transform is similar to the use of PDE for image
analysis, which is an important ﬁeld. Nonlinear PDEs play a crucial role there. What has
been demonstrated in the present work is that when PDEs are used for conservation law
problems, appropriate high order PDEs are indispensable. What has not been examined
in the present work is the role of nonlinear geometric PDEs for conservation law problems.
This aspect is under our consideration. We hope that our work will lead to new exploration
in the use of PDEs for systems of hyperbolic conservation laws.

90

Chapter 3
High order fractional PDE transform
for molecular surface construction

3.1

Motivations

Recently, great interests have been witnessed in fractional calculus modeling in many ﬁelds
of sciences and engineering, from geophysics to biology. Fractional derivatives extend the
concept of ordinary derivatives and serves as a good tool for taking into account memory
mechanism in the random walk and anomalous diﬀusion in physical problems [125, 36]. In
fact, a fractional derivative occurs as common as an exponent. For example, it is well known
that the Kolmogorov scaling exponent that predicts turbulent energy spectrum follows a
power law E(w) ∝ w−p , where w is the wavenumber, p is a real number, and 0 < p < 3.
Such a power law in the wavenumber corresponds to the fractional derivative of order α
in the coordinate space, since the inverse Fourier transform of an exponent gives rise to a
derivative of the same order. Fractional derivative has found its success in physical and
biological modeling [8, 197, 22], ﬁnancial analysis [74, 136, 142], and image processing [9].
Currently, most attention in the ﬁeld is paid to the fractional derivatives of order less than
2. high-order fractional derivatives are hardly used, partly due to the limited understanding
of their physical meanings.
In our work, we introduce arbitrarily high-order fractional PDEs and the fractional PDE
91

transform for the analysis of molecular surfaces in molecular biosciences. Since arbitrarily
high integer order nonlinear PDEs was originally introduced to account for hyperdiﬀusion in
physical and biological systems [178], it is natural to consider arbitrarily high-order fractional
PDEs and the fractional PDE transform, which utilizes arbitrarily high-order fractional
derivatives. It is well known that fractional derivatives can be deﬁned in a few diﬀerent
ways. We discuss the ﬁnite diﬀerence, the Fourier representation and the integral forms of
fractional derivatives. A variational derivation of fractional PDEs is presented. To avoid the
strict stability constraints of solving high-order fractional PDEs and to achieve the desirable
eﬃciency in the surface generation, we make use of the fast fractional Fourier transform
(FFFT) algorithm to realize the fractional PDE transform. The present algorithm is of
order O(N ln N ). Extensive numerical test and application validate the proposed fractional
PDE transform.

3.2

Theory and algorithm

This section provides theory and algorithm for arbitrarily high-order fractional PDEs and
fractional PDE transform. To establish notations, we ﬁrst brieﬂy discuss a few approaches
for fractional derivatives. The fractional PDE transform is constructed via a variational
analysis.

3.2.1

Fractional derivatives

There are many diﬀerent approaches for the deﬁnition and understanding of fractional derivatives. Among them, fractional ﬁnite diﬀerence schemes, fractional Fourier representations
and the integral forms of fractional derivatives are the most popular ones.

92

3.2.2

Fractional ﬁnite diﬀerence schemes

Finite diﬀerences provide a natural representation of fractional derivatives. There are three
types of standard ﬁnite diﬀerence schemes, namely, forward, backward and central ﬁnite
diﬀerences. As a consequence, three fractional ﬁnite diﬀerence schemes can be deﬁned accordingly. For any integer α > 0, the α-th order backward ﬁnite diﬀerence operator is given
by [125]


α

∆α u(x) =
h,b



 α 
(−1)m 
 u(x − mh),
m=0
m

h > 0,

(3.1)

where the Binomial coeﬃcients




α!
 α 
.

=
m!(α − m)!
m

(3.2)

For example, when α = 1 and 2, the ﬁrst and second order backward ﬁnite diﬀerence schemes
are

∆h,b u(x)
u(x) − u(x − h)
=
h
h

(3.3)

and
∆2 u(x)
h,b
h2

=

u(x) − 2u(x − h) + u(x − 2h)
,
h2

(3.4)

respectively.
Equation (3.1) can be generalized to the case for any real number α > 0, in which the

93

Binomial coeﬃcients are expressed in terms of Gamma function Γ, i.e.,




Γ(α + 1)
 α 
,

=
Γ(m + 1)Γ(α − m + 1)
m
and hence
a2 −a1
h

∆α u(x) =
h,b

(−1)m
m=0

Γ(α + 1)
u(x − mh),
Γ(m + 1)Γ(α − m + 1)

where a1 and a2 are the lower and upper bound of diﬀerentiation, respectively. Thus, the
factional derivative of order α can be deﬁned as

u(α)

= lim

∆α u
h,b

h→0

hα

.

(3.5)

One special property of the fractional ﬁnite diﬀerence operator in Eq. (3.1) is that it utilizes
all the node points to the left. Similarly, one can deﬁne the fractional derivatives of α-th
order by forward and central ﬁnite diﬀerence operators. Finite diﬀerence based fractional
derivatives ﬁnd their applications in anomalous diﬀusion [126].

3.2.2.1

Fractional Fourier schemes

It is well known that Fourier transform provides a practical approach for computing integer derivatives and solving diﬀerential equations. It can be similarly used for fractional
derivatives [125]. Consider the Fourier transform
∞

u(w) =

e−iwx u(x)dx.

−∞

94

(3.6)

The Fourier transform of the backward ﬁnite diﬀerence operator in Eq. (3.1) is given by


∞

∆α u(w) =
h,b



 α  −imhw
(−1)m 
u(w) = (1 − e−ihw )α u(w).
e
m=0
m

(3.7)

where Binomial coeﬃcients vanish when m > α. Consequently, the corresponding Fourier
transform for (3.5) is

u(α) (w)

= lim

∆α u
h,b

h→0 hα

= lim

h→0

1 − e−ihw
h

α

u(w) = (iw)α u(w),

(3.8)

where the last step makes use of the Taylor expansion e−ihw = 1 − ihw + O(h2 ) as h → 0.
There are many other ways to derive Eq. (3.8). Consider the central ﬁnite diﬀerence in
the Fourier domain deﬁned as [9]

Dn u = (eihw/2 − e−ihw/2 )Dn−1 u = · · · = (eihw/2 − e−ihw/2 )n u,

(3.9)

where the nth-order ﬁnite diﬀerence Dn u is deﬁned by

Dn u(x) = D(Dn−1 u) = · · · = Dn−1 (Du) = Dn−1 u x +

h
2

−u x−

h
2

. (3.10)

As an extension of the Fourier transform of the nth-order ﬁnite diﬀerence, one has

Dα u = eihw/2 − e−ihw/2

95

α

u.

(3.11)

Consequently, generalized fractional-order derivative is given by

u(α)

Dα u
= lim
= lim
h→0
h→0 h

eihw/2 − e−ihw/2
h

α

u = (iw)α u,

(3.12)

where the last step uses Taylor expansion eihw/2 = 1 + ihw/2 + O(h2 ) and e−ihw/2 =
1 − ihw/2 + O(h2 ) when h → 0. As a result, Eq. (3.8) is also concluded.
When α is an integer, i.e., α = n, Eq. (3.8) becomes

u(n) = (iw)n u,

(3.13)

which is consistent with the classical derivative in the Fourier representation. Moreover,
when α = 2m is an even number, we have

u(2m) = (−1)m (w)2m u.

(3.14)

Furthermore, using the Euler’s formula i = exp(iπ/2), Eq. (3.8) can be reformulated as

u(α) = eiαπ/2 wα u.

(3.15)

The expression in Eq. (3.8) can be numerically computed via Eq. (3.15).

3.2.2.2

Integral forms

For any real number α > 0, the Caputo fractional derivative of order α is deﬁned as [29]

Dα u =

x
1
u(n) (t)
dt,
Γ(n − α) 0 (x − t)α+1−n

96

n − 1 < α < n,

(3.16)

where it requires the absolute integrability of u up to the n-th integer derivatives. Meerschaert has developed a discrete approximation to integral form of fractional derivative in
Eq. (3.16) based on the shifted Gr¨nwald approximation [126].
u
Somewhat less restrictive integral deﬁnitions have also been given in the literature by
making use of Mittag-Leﬄer-type functions [121]. For any real number β > 0, denote by
I β f (x) the Riemann-Liouville integral
t2
1
β
It ,t u(x) =
u(t)(x − t)β−1 dt.
1 2
Γ(β) t1

(3.17)

Then the left and right Riemann-Liouville fractional derivatives of order α for a function
u(x) are given by [2]

α
Dl u =

dn n−α
I
u(x)
dxn a,x

(3.18)

and

α
Dr u = (−1)n

dn n−α
I
u(x),
dxn x,b

(3.19)

respectively. In the integral, a and b are the lower and upper bounds, which are typically
taken as −∞ and ∞, respectively. The right Riemann-Liouville fractional derivative operates
on the right (or future) state of the function u(x) and thus is less popular in the literature.
However, it is useful in variational formulations.
If α is an integer, these derivatives are deﬁned in the usual sense, namely

α
Dl u =

dα
u(x),
dxα

α = 1, 2, · · · ,
97

(3.20)

and

α
Dr u = (−1)α

dα
u(x),
dxα

α = 1, 2, · · · .

(3.21)

The left Riemann-Liouville fractional derivative is one of the most popular deﬁnitions. Some
interesting properties and history of Riemann-Liouville fractional derivatives have been discussed in Ref [121].

3.2.3

High-order fractional PDEs and fractional PDE transform

3.2.3.1

Variational derivation of high-order fractional PDEs

Motivated by the PDE transform and fractional derivatives, we propose a fractional PDE
transform. A key component of the fractional PDE transform is an arbitrarily high-order
fractional PDE, which is deﬁned via the fractional variational principle. Many authors
have discussed the fractional variational principle [2, 9]. For any α ∈ R and 0 < α < ∞,
denote by

α

=

∂lα ∂lα ∂lα
,
,
∂xα ∂y α ∂z α

a gradient vector in R3 , where

∂lα
is the left Riemann∂xα

Liouville fractional derivative of order α in x, according to Eq. (3.18); similar notations
are adopted for the fractional derivatives in y and z. The adjoint of
α
∗

=

α

is deﬁned as

α
α
α
∂r ∂r ∂r
, α , α , with the right Riemann-Liouville fractional derivative in (3.19).
∂xα ∂y ∂z

We now consider the general energy functional

E(u,

α u,

α u)
∗

=

F (u,

α u,

α u) dr
∗

=

Λ | α u|2 , | α u|2 + (X − u)2 dr,
∗
(3.22)

where Λ and (X − u)2 are similar to the corresponding terms in Eq. (2.7). The Euler-

98

Lagrange equation of (3.22) is [2]:

∂F
+
∂u

α
∗

·

∂F
∂

αu

+

α

·

∂F

= 0.
∂ αu
∗

(3.23)

Applying (3.23) to (3.22), we obtain

1 α l
1 α r
· Λuα + (X − u) = 0,
∗ · Λuα +
2
2

(3.24)

where
Λl = −∂Λ/∂ α u
uα

and

Λr = −∂Λ/∂ α u.
uα
∗

(3.25)

Note that Eq. (3.24) is in the similar fashion as Eq. (2.11).

Remarks
1. We can introduce an artiﬁcial time t and convert Eq. (3.24) into a time-dependent
fractional PDE
1 α l
1 α r
∂u
=
· Λuα + (X − u),
∗ · Λuα +
∂t
2
2

(3.26)

then the fractional PDE transform is constructed by Eq. (2.14) with the solution operator
L deﬁned via Eq. (3.26).
2. Furthermore, setting Λ = Λ(µ, ν) with µ = | α u|2 , ν = | α u|2 , we can apply the
∗
chain rule to (3.25), i.e.,

Λl = −
uα

∂Λ ∂µ
∂Λ α .
= −2
u = −2Λµ α u
αu
∂µ ∂
∂µ

99

(3.27)

and
Λr = −
uα

∂Λ ∂ν
∂Λ α .
α
α u = −2 ∂ν ∗ u = −2Λν ∗ u.
∂ν ∂ ∗

(3.28)

Substituting (3.27)-(3.28) into (3.26), we obtain

∂u
= − α · Λµ α u −
∗
∂t

α

· (Λν α u) + (X − u),
∗

(3.29)

where Λµ | α u|2 , | α u|2 , Λν | α u|2 , | α u|2 can be viewed as anisotropic diﬀusion co∗
∗
eﬃcients.
3. A special case is that when only

αu

(or

α u)
∗

appears in (3.22) and Λ takes the

simple linear function Λ(µ) = µ, Eq. (3.29) turns into

∂u
=− α·
∗
∂t

αu +

(X − u).

By (3.20) and (3.21), we can recover the simple heat equation

∂u
=
∂t

2u +

(X − u) for

α = 1.

3.2.3.2

Fractional hyperdiﬀusions derivation of high-order fractional PDEs

In this section, we derive from Eq. (3.29) an arbitrarily high-order fractional hyperdiﬀusion
equation which has the similar structure of the hyperdiﬀusion equation (2.3). Actually,
letting
m

Λ = Λ(µ),

| αq u|2 ,

µ=
q=1

100

where m ∈ N, αq ∈ R, and 0 < α1 < α2 < ... < αm < ∞, we have, according to Eq. (3.29),

m
 ∂u



=−
 ∂t



q=1


αq
∗



m

· jαq + e u,

| αq u|, t ,
q=1

(3.30)







 αq
 j = d (u, | αq u|, t) αq u,
αq


m





| αq u|, t = (X − u) and dαq u,

where e u,
q=1



m



| αq u|, t = Λµ 
q=1



m

| αq u|2 .

q=1

Setting m = 1, αm = 1, e(u(r, t), | u(r, t)|, t) = 0, and recalling (3.20)-(3.21), one can
reproduce the well-known Perona-Malik equation for anisotropic diﬀusion [134]:

∂u
=
∂t

· [dα (| u|, t) u] .

In the current work we study the high-order fractional hyperdiﬀusion equation with the
linear form of Λ, i.e, Λ(µ) = µ, and neglect the term e(u, | αq u|, t).

3.2.4

Numerical algorithms

In our previous work, high-order nonlinear evolution PDEs have made vital contributions
to mode decomposition [171], image analysis [169], and molecular surface generation [204].
However, solving high-order PDE in Eq. (3.26) directly in physical domain can be numerically diﬃcult for some practical applications. A major problem in solving high-order nonlinear
evolution PDEs is their stability constraints: in general, in explicit numerical methods, the
time increment ∆t needs to be proportional to the n-th power of the spatial discretization
h, i.e., ∆t ∼ hn , for an n-th order PDE. Although the constraint has little impact to most
image analysis, in which the spacing h is unit, it does lead to some diﬃculties in other appli101

cations. Typically, alternating direction implicit (ADI) methods are implemented to by-pass
the stability constraint in solving high-order PDEs [14]. Another approach used in previous
PDE transform is the Fourier pseudospectral method via the fast Fourier transform (FFT).
These techniques can be analogically applied to the fractional high-order PDE transform. In
such an approach, one has to linearize Eq. (3.30) and then solves original nonlinear equation
by an iterative procedure as shown in the part of PDE transform.
We consider a linearized fractional PDE of (3.30), i.e. Λ(µ) = µ, µ = | α u|2 , and m = 1,
then

∂u
= −dq α ·
∗
∂t

α u,

α ∈ R+ ,

(3.31)

where dq is a coeﬃcient independent of u. Denote by F[u](iw) the discrete Fourier transform
of u(r, t), t ≥ 0.
When Eq. (3.31) is subject to initial value u(r, 0) = X k and periodic boundary conditions, it can be easily solved in the Fourier representation. The Fourier transform of Eq.
(3.31) is given by

∂F[u](w)
= −dq F[ α ·
∗
∂t

α u](w),

(3.32)

Because of relations

α
F[Dl u] = (iw)α F[u](w),

α
F[Dr u] = (−iw)α F[u](w),

102

(3.33)

Eq. (3.32) turns into

∂F[u](w)
= (−1)α+1 i2α (w)2α dq F[u](w) t ≥ 0,
∂t

(3.34)

where i2α is evaulated as i2α = eiαπ . Therefore, we iteratively solve the linearized high-order
fractional PDE until a predesigned accuracy is reached. We can express the solution of the
original nonlinear high-order PDE as [169]

k

ˇ
X (w) = LX k (w)

(3.35)

ˇ
ˇ
where X and X k are the Fourier transforms of X k and X k , respectively. Here Fourier transform of the evolution operator is L = exp (−1)α+1 i2α (w)2α dq t . The solution algorithm
developed in the above Fourier domain is called a fast PDE transform.
Given initial value u(r, 0) = X k , the evolution PDE (3.31) is solved by the fast Fourier
transform (FFT) as shown in Eq. (3.35). Speciﬁcally, after computing u = FFT(u) by the
FFT, we apply the fractional PDE transform in the Fourier space given by Eq. (3.35) so as
to obtain u = Lu with corresponding propagation time t. Finally, the inverse FFT (IFFT) is
ˇ
implemented for u to produce u = IFFT(u). As a result, the new value of u after propagation
ˇ
ˇ
time t, i.e., u(r, t), will be given by u. In our experiments, the FFT is implemented by the
‘ﬀtw’ software package (www.ﬀtw.org).
It is noted that since the fractional PDE transform is implemented in Fourier space, the
term (iw)α will involve the complex number when α is not an even integer. Moreover, the
complex u will result in complex u. However, the information of the biomolecular surface is
ˇ
preferably represented by the real number as the isosurface is from the value of u(r, t), which

103

should be real. In the practical implementation, one can either take the real component of
Eq. (3.34) (see Ref. [9] for a similar treatment) or the amplitude of the term (iw)α . In this
work, we explore the performance of both approaches.

3.3

Numerical test and validation

In this section, the proposed fractional PDE transform is tested and validated by molecular
surface generation. The initial condition of Eq. (3.31) is set as initial biomolecular surface
data and the isosurface obtained from the solution u(r, t) at time t will be a desired molecular
surface after the PDE transform.

3.3.1

Initial data for evolution PDEs

The initial value of Eq. (3.31) is set to be a molecular hypersurface function, which contains
atomic coordinates and radius of a molecule. In our computation, three types of initial data
are used.

3.3.1.1

Type I: Piecewise-constant initial values

Piecewise-constant initial values are commonly used in our earlier work for molecular surface
generation [15, 16, 37]. In this work, we set the piecewise-constant initial value as

u(r, 0) =



 0,

r∈


 1,

i O(ri , ri );

(3.36)

otherwise,

where ri and ri are the atomic coordinate and radius of atom i, respectively. The union
of all atomic spheres O(ri , ri ) is denoted by
104

i O(ri , ri ).

This type of piecewise-constant

initial value is diﬀerent from the initial value in our geometry ﬂow based solvation model
[37], where the two constants 0 and 1 are switched.

3.3.1.2

Type II: Maximum Gaussian initial values

The second type of initial values is deﬁned by the maximum of Gaussian functions [68, 204]
as follows

u(r, 0) = maxi s exp −

2
r − ri 2 −ri
2
re

.

(3.37)

In this work, we set s = 1 and re = 3 as recommended in the literature [68].

3.3.1.3

Type III: Summation Gaussian initial values

The third type of initial values is taken from the summation of Gaussian functions [21, 80, 64],

2
exp −κ( r − ri 2 −ri ) .

u(r, 0) =

(3.38)

i

In our simulation, we set the decay rate κ = 1 [198].
The initial molecular surface can be obtained as the isosurface by setting u(r, 0) = S,
where S is a pre-designated isovalue of u, for example, u(r, 0) = 1 [198, 204]. Although the
two types of Gaussian initial data are able to oﬀer smoother isosurfaces, they may still result
in isosurfaces with singularities as demonstrated in the following sections.

3.3.2

Eﬀects of fractional order and propagation time

It is noted from the evolution PDE in Eq. (3.31) that, the performance of the proposed
fractional PDE transform depends on the fractional order α, propagation time t, and diﬀusion
105

coeﬃcient dq . In the present work, we mainly focus on the role of fractional order and
propagation time while set dq = 1.
Figure 3.1 shows the Fourier response of the fractional PDE transform, which is equivalent
to the operator L in Eq. (3.35) at diﬀerent orders. The overall amplitudes of the Fourier
responses of fractional PDE transforms with diﬀerent orders are given in Figure 3.1(a),
and 3.1(b)-(c) show the corresponding real and imaginary components, respectively. It can
be concluded that the amplitude of Fourier response behaves like classic ﬁlters. However,
the real component has unfavorable negative values, and the magnitude of the imaginary
component is relatively small. These results provide an understanding of the behavior of
high-order fractional PDEs.
Figure 3.1(a) also depicts the eﬀects of fractional PDE transform with diﬀerent orders
(α = 1.5, 5.5, 11.5). With a ﬁxed propagation time t = 102 , the fractional PDE transform
of order α = 1.5 oﬀers a low-pass ﬁlter for very low frequency components. Whereas,
the fractional PDE transform of order α = 11.5 is able to pass a much larger range of
frequency components. Based on this feature we choose diﬀerent fractional orders in mode
decompositions and other applications of high-order PDEs.
The eﬀect of propagation time, which is equivalent to the impact of the nonlinear coeﬃcient in the PDE transform, is shown in Figure 3.2, with ﬁxed fractional orders α ∈ (11, 12)
and diﬀerent propagation time. From Figure 3.2 (a) to Figure 3.2 (c), there display the real
components of the Fourier response of the fractional PDE transform with order α = 11.2,
α = 11.5 and α = 12.0, respectively. It is concluded that the longer the propagation time is,
the faster is the damping of real component. Thus we can conjecture that in order to obtain
a smoother molecular surface, a longer propagation time is preferred. We also need to keep
in mind that the molecular surface can not be over-smoothed or over-smeared, otherwise it
106

α=1.5
α=5.5
α=11.5

Amplitude

1
0.8
0.6
0.4
0.2
0
0

1
2
Wavenumber
(a)

α=1.5
α=5.5
α=11.5

1
0.8
Real

3

0.6
0.4
0.2
0
0

1
2
Wavenumber
(b)

α=1.5
α=5.5
α=11.5

1
Imaginary

3

0.8
0.6
0.4
0.2
0
0

1
2
Wavenumber
(c)

3

Figure 3.1: Frequency response of Eq. (3.34) with t = 102 and diﬀerent orders. (a) Over all
amplitude; (b) Real component; (c) Imaginary component. Red: α = 1.5; blue: α = 5, 5;
and purple: α = 11.5.

107

will become unphysical. It is observed that when the fractional order is closer to an odd
number, the real components become more and more oscillatory with all propagation time.
Furthermore, Figure 3.3(a) shows the visualization of real component of Fourier response in
2D for α = 11.5 and t = 102 ; Figure 3.3(b) depicts its contours. Since there are negative
oscillations in the Fourier response, there are corresponding oscillations in contours as shown
in Figure 3.3(b).

3.3.2.1

Test on a three-atom system

We start the test cases with a simple three-atom system, in which the initial atomic coordinates are given as (0, 0, 1.8), (0, 0, −1.8), (0, 3.12, 0) with a uniform radius of 1.8. Figure 3.4
shows the isosurface of three-atom system from piecewise constant initial value (Figure 3.4
(a)) and after fractional PDE transform by diﬀerent fractional orders α = 1.5 (Figure 3.4
(b)), 5.5 (Figure 3.4(c)) and 11.5 (Figure 3.4(d)); the propagation time is t = 102 and the
isovalue is taken as u(r, t) = 1. We conclude that the piecewise constant initial value does
not give an acceptable molecular surface and the low order transform (α = 1.5) over-smears
the molecular surface. Similar problem of a low order PDE transform at α = 2 was also
found by Zheng et al [204]. However, both high-orders transforms (α = 5.5 and α = 11.5)
render satisfying molecular surfaces. In the work of integer PDE transform based biomolecular surface generation [204], although α = 6.0 works well for small molecular systems, it
was found that α = 6.0 over-smears large biomolecules such as proteins. In the following
test cases, we use the fractional order α = 11.5.

108

Real

1

6

t=10
4
t=10
t=102

0.5
0
−0.5
0

1
2
Wavenumber
(c)

6

1

t=10
t=104
t=102

0.8
Real

3

0.6
0.4
0.2
0
0

1
2
Wavenumber
(b)

t=106
4
t=10
2
t=10

1
0.8
Real

3

0.6
0.4
0.2
0
0

1
2
Wavenumber
(c)

3

Figure 3.2: Real component of frequency response of Eq. (3.34) in 1D, with diﬀerent values
of α and t: (a) α = 11.2; (b) α = 11.5; (b) α = 12.0. Red: t = 106 ; blue: t = 104 ; and
purple: t = 102 .

109

1

0.5

0

−0.5
−2

0

2

2

0 −2

(a)
3
2

0.8

1

0.6

0
0.4
−1
0.2
−2
0

−3
−2

0

2

(b)
Figure 3.3: Real component of Fourier response of Eq. (3.34) in 2D for α = 11.5 and t = 102 .
(a) function value view; (b) contour view.

110

(a)

(b)

(c)

(d)

Figure 3.4: Isosurfaces of a three-atom system after the fractional PDE transform of diﬀerent
fractional orders with the piecewise constant initial value. The isovalue is taken as u(r, t) = 1
and propagation time is t = 102 . (a) Initial value; (b) α = 1.5; (c) α = 5.5; (d) α = 11.5.
3.3.2.2

Test on a four-atom system

We next test the eﬀects of diﬀerent propagation time and choices of initial data through a
four-atom system, in which the initial atomic coordinates are given as (0, 0, 1.8), (0, 0, −1.8),
(0, 3.12, 0) and (0, −3.12, 0), with a uniform radius of 1.8. The fractional PDE transforms
of diﬀerent propagation time with three types of initial data discussed in section 3.3.1, are
illustrated in Figures 3.5, 3.6 and 3.7, respectively, with the isovalue u(r, 0) = 1 and order
α = 11.5. For the initial molecular surfaces, the one from Type I is non-smooth (Figure
3.5(a)); isosurface of Type II has singularities (Figure 3.6(a)); and Type III is smooth (Figure
3.7(a)). However, it should be noted that the three types of initial values have diﬀerent
ranges: initial data of Type I is restricted within the range [0, 1], while the values of Types II
and III can be greater than one. Although the isosurface of Type III is smooth at u(r, 0) = 1,
the surface from other isovalues may not be. For example, the surface obtained by setting
u(r, 0) = 2, as showed in Figure 3.7(d), has obvious singularities.
It is noted that Type II initial value depends on the parameter re and Type III data
relies on the decay rate κ. Figure 3.8 shows the isosurfaces from Type III initial data at
two diﬀerent decay rates, κ = 1 and κ = 1/9. It is observed that κ = 1 gives a much

111

(a)

(b)

(c)

Figure 3.5: Isosurface of a four-atom system after fractional PDE transform of order α = 11.5
and diﬀerent propagation time, with the piecewise constant initial value. The isovalue is
taken as u(r, t) = 1 (a) initial value; (b) t = 1; (c) t = 104 .
better resolution of the molecular surface than κ = 1/9 does; this fact is consistent with the
previous results [198].
The results after fractional PDE transform with order α = 11.5 after diﬀerent propagation
time are illustrated in the rest subﬁgures of Figures 3.5, 3.6 and 3.7. Although the piecewise
constant (Type I) initial data is non-smooth, it can oﬀer smooth surfaces by the fractional
PDE transform after certain proper propagation time, as noted in Figures 3.6(b) and (c).
On the other hand, the singularities resulting from Type II initial data can also be removed
at a long propagation time t = 104 but not after a short time t = 1. Moreover, it is diﬃcult
to ﬁlter the singularities in the isosurface u(r, t) = 2 with Type III initial value. One may
also pay attention to the time cost for generating three types of initial surfaces. Since Types
II and III involving searching the maximum or taking the summation, they require more
CPU time than Type I does. Indeed, for the four-atom system at h = 1, the computational
results show that the CPU cost for surface generation using Type I initial data is about
0.001 seconds, while about 0.2 seconds for using Types II and III initial data.

112

(a)

(b)

(c)

Figure 3.6: Isosurfaces of a four-atom system after fractional PDE transform of order α =
11.5 and diﬀerent propagation time, with maximum Gaussian initial value. The isovalue is
taken as u(r, t) = 1. (a) Initial value; (b) t = 1; (c) t = 104 .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3.7: Isosurfaces of a four-atom system after fractional PDE transform of oder α = 11.5
and diﬀerent propagation time, with the summation Gaussian initial value. The isovalue is
taken as u(r, t) = 1 for (a)-(c) and u(r, t) = 2 for (d)-(f). (a) Initiale value; (b) t = 1; (c)
t = 104 ; (d) initial value; (e) t = 1; (f) t = 104 .

113

(a)

(b)

Figure 3.8: Isosurface of a four-atom system with initial value of Type III (Summation
Gaussian initial values) at diﬀerent κ. The isovalue is taken as u(r, 0) = 1. (a) κ = 1; (b)
κ = 1/9.
3.3.2.3

Test on two proteins

The fractional PDE transform method is further validated on realistic proteins. In our
numerical experiments, the protein data, i.e., atomic coordinates and radii, are obtained
from the PDB bank (www.pdb.org) and prepared with the PDB2PQR package [50]. Two
proteins, PBD IDs 1R69 and 1FCA, are used in our test. The results are depicted in Figures
3.9 and 3.10, respectively. Initial values of Type I are used for both proteins (Figure 3.9(a)
and 3.10(a)). Diﬀerent orders (α = 1.5, 11.5) and the same propagation time (t = 104 )
are applied to the protein 1R69, see Figure 3.9(b)-(c), while diﬀerent propagation time
(t = 1, 104 ) and the same order (α = 11.5) are for the protein 1FCA, see Figure 3.10(b)-(c).
It can be concluded that the transform with relatively longer propagation time and higher
order produces better molecular surfaces.
Based on the tests on the model three/four-atom systems and the realistic proteins, we
summarize that even with the simplest Type I initial surface, satisfactory molecular surfaces
can be obtained if the parameters are chosen properly. Speciﬁcally, we use order α = 11.5,
propagation time t = 104 and isosurface u(r, t) = 1 for later applications unless other values
114

are indicated. The initial values are all of Type one for connivence.

(a)

(b)

(c)

Figure 3.9: Surface of a protein (PDB code: 1R69) after fractional PDE transform of propagation time t = 104 and diﬀerent orders. The isovalue is taken as u(r, t) = 1. (a) Initial
value (Type I); (b) α = 1.5; (c) α = 11.5;

(a)

(b)

(c)

Figure 3.10: Surface of a protein (PDB code: 1FCA) after fractional PDE transform of order
α = 11.5 at diﬀerent propagation time. The isovalue is taken as u(r, t) = 1. (a) initial value
(Type I); (b) t = 100 ; (c) t = 104 ;

3.3.3

Computational eﬃciency

exit
Computational eﬃciency is a crucial issue in molecular surface generation. Surface generation time is a bottleneck for many practical applications such as the Poisson-Boltzmann
115

based molecular dynamics [66] and the solution of coupled nonlinear Poisson-Nernst-Planck
equations [202, 201]. Here we examine the computational eﬃciency of the fractional PDE
transform based biomolecular surface generation. Table 3.1 compares the CPU costs (in
seconds) for surface generation in Cartesian meshes between the fractional PDE transform
(FPDE) with α = 11.5, t = 104 and the MSMS approach [143] at a typical triangle density
10. In the table, Na is the total number of atoms in the protein after the PDB2PQR conversion. The total time Ttot of the MSMS approach includes the CPU cost from generating
surface in the 2D Lagrangian representation Ts and the time cost of converting it to the
3D Eulerian representation Tc , which is required in many applications [204]. Essentially, the
CPU time used to create a triangular representation of a protein surface, i.e., Ts , is similar to
the total CPU time for the the present fractional PDE transform to generate a 3D Eulerian
representation of the protein surface. However, the additional CPU used in the Eulerian conversion is quite large, which causes problems in our earlier molecular dynamics simulations
[66]. Moreover, the MSMS method encounters problem in generating the triangular surface
meshes for three relatively large proteins, 2NCD, 1IWO and 2E2J at the designed density
10. The present fractional PDE transform is very stable and robust for large proteins.
It is interesting to note that the computational cost scales essentially linearly with respected to the number of atoms in proteins for the present fractional PDE transform approach.

3.3.4

Surface areas and surface enclosed volumes

To quantitatively characterize the surface generated by the fractional PDE transform, we
examine the surface areas and their enclosed volumes for a number of proteins. These
quantities are frequently used in biophysical modeling, such as the solvation free energy
estimation.
116

PDB-ID
1R69
1A2S
1SH1
2PDE
1VII
1A7M
1VJW
1FCA
1HPT
1MBG
1A63
1SVR
2NCD
1IWO
2E2J

Na
596
667
702
729
828
858
903
997
1272
1435
2065
2809
5665
30870
58046

Ts
by MSMS
0.49
0.43
0.56
0.46
0.54
0.64
0.65
0.62
0.92
0.96
1.38
1.58
N/A
N/A
N/A

h=1
MSMS
FPDE
Tc
Ttot
Ttot
2.84 3.33
0.05
3.00 3.43
0.06
3.13 3.69
0.09
2.91 3.37
0.11
3.03 3.57
0.06
3.73 4.37
0.05
3.59 4.24
0.08
3.52 4.14
0.08
5.11 6.03
0.07
5.15 6.11
0.11
7.92 9.30
0.19
8.73 10.31 0.19
N/A N/A
0.94
N/A N/A
3.20
N/A N/A
3.48

h = 0.5
MSMS
FPDE
Tc
Ttot
Ttot
2.93 3.42
0.32
3.25 3.68
0.53
3.27 3.83
0.43
3.31 3.77
1.28
3.82 4.36
0.60
3.90 4.54
0.57
3.56 4.21
0.67
3.72 4.34
0.55
5.21 6.13
1.08
5.34 6.30
1.21
8.09 9.47
1.20
9.29 10.87 1.82
N/A N/A
9.35
N/A N/A 23.97
N/A N/A 31.36

Table 3.1: Comparison of CPU costs (in seconds) for surface generation in Cartesian meshes
by the fractional PDE transform (FPDE) with α = 11.5, t = 104 and the MSMS approach.
Na : total number of atoms; Ts : CPU time for MSMS to generate surface in 2D Lagrangian
representation; Tc : CPU time for converting the Lagrangian representation to 3D Eulerian
representation. Ttot : Total CPU cost for FPDE and MSMS; h: grid resolution.
The results of fractional PDE transform are in the Eulerian representation, i.e., a 2D
surface of a biomolecule embedded in a 3D Cartesian grid. One has to extract the surface
area and the enclosed volume by appropriate computational algorithms. A second-order
surface integration scheme has been developed in our earlier work [66]

f (x, y, z)dS ≈
Γ

f (xo , yj , zk )
(i,j,k)∈I

| ny |
| nz |
| nx |
+ f (xi , yo , zk )
+ f (xi , yj , zo )
h
h
h

h3
(3.39)

where h is the grid resolution, (xo , yj , zk ) is the intersecting point of the surface Γ and
the meshline that passes through the point (xi , yj , zk ), and nx is the x component of the
unit normal vector at (xo , yj , zk ). Similar notations are used for y and z directions. The
117

PDB-ID
1VJW
1FCA
1SH1
1R69
1MBG
1A2S
1SVR
1VII
2PDE
1A7M
1A63
1HPT

Area
2847.6
2723.2
2680.7
2978.9
2980.2
4316.1
4582.6
2402.0
2659.2
7514.9
6884.9
3182.9

FPDE
Volume
7626.1
6857.1
6418.2
8817.1
7896.8
11629.6
11990.2
5091.7
6080.3
24125.5
18168.8
7706.9

∆Gp
-1236.0
-1206.0
-746.9
-1081.9
-1337.9
-1902.1
-1699.0
-888.3
-805.3
-2129.9
-2347.0
-807.5

Area
2785.8
2546.9
2743.8
3054.5
3070.7
4436.8
4644.0
2476.2
2715.9
7733.6
6973.7
3262.9

MSMS
Volume
7693.2
7028.1
6400.7
8777.2
7843.3
11569.7
11961.9
5050.7
5994.6
24025.1
18121.1
7660.9

∆Gp
-1237.9
-1200.1
-753.3
-1089.5
-1346.1
-1913.5
-1711.2
-901.5
-820.9
-2155.5
-2373.5
-811.6

Table 3.2: Comparison of surface areas (unit: ˚2 ), surface enclosed volumes (unit: ˚3 ) and
A
A
electrostatic free energies of solvation ∆Gp (unit: kcal/mol) with surfaces generated by the
fractional PDE transform (FPDE) and the MSMS approach. The grid resolution is 0.5 ˚ .
A
summation involves irregular grid points which are the points with at least one neighbor
from the other side of the surface Γ in the second-order ﬁnite diﬀerence discretization. Here
I is the set of irregular grid points inside or on the interface [66]. Formula (3.39) is used for
surface area calculation by setting f = 1.
The volume integral of a function f is approximated by



f (x, y, z)dr ≈
Ωm

1
2

f (xi , yj , zk )h3 +

f (xi , yj , zk )h3 

(3.40)

(i,j,k)∈J2

(i,j,k)∈J1

where J1 is the set of the points inside Ωm and J2 = J1

JIrr , where JIrr indicates the set

of irregular points. The volume is obtained by setting f = 1 in Eq. (3.40).
Table 3.2 lists surface areas and surface enclosed volumes of protein surfaces generated
by both the fractional PDE transform and the MSMS approach. In this test, the fractional
PDE transform is of order 11.5, propagation time t = 104 , and the grid resolution is 0.5 ˚.
A
It can be concluded that the present results are quite consistent with those obtained from
118

the MSMS package, thus the proposed PDE transform is useful in practical applications.

3.4

Applications

Solvent-solute interface plays a critical role in implicit solvent models for applications to
electrostatic calculation, solvation analysis and molecular dynamics simulations. It is also
important in variational multiscale models for biomolecular systems [180, 184]. In this section, we ﬁrst apply the fractional PDE transform based biomolecular surfaces to electrostatic
calculation of a group of proteins and then analyze the solvation free energies of a set of small
compounds, for which have the experimental data are available.

3.4.1

Electrostatic analysis of 15 proteins

3.4.1.1

Electrostatic potential

Once the surface has been generated, the electrostatic potential can be evaluated according
to the Poisson equation [194, 67, 33]

− ·(ε φ) =

qj δ(r − rj ),

(3.41)

j

where qj are the (fractional) charges of atoms at position rj (j = 1, 2, ..., Na ). The dielectric
function

is deﬁned as a piecewise constant function

ε(r) =



ε , r ∈ Ω ,
 m
m


 εs ,

119

r ∈ Ωs ,

(3.42)

where εm = 1 and εs = 80 are the dielectric constants in the molecular and solvent regions,
respectively; these two regions are separated by the molecular surface. Equation (3.41)
is a typical elliptic interface problem with discontinuous coeﬃcients and singular sources.
It is very diﬃcult to construct second-order convergent methods for this equation in the
biomolecular context due to the geometric complexity, complex interface, singular charge
sources and geometric singularities [194, 67, 33]. In this work, we make use of the secondorder convergent matched interface and boundary (MIB) method [196, 195, 199, 208, 207, 33]
to solve Eq. (3.41).
Figure 3.11 displays the electrostatic surface potential of three proteins obtained from
solving the Poisson equation (3.41). These potential values are projected on the protein
surfaces generated by using the proposed fractional PDE transform. In the ﬁgures, red
and blue colors represent negative and positive surface electrostatic potentials, respectively.
Surface electrostatic potential is highly correlated with ligand binding domains, which are
important to drug design.

(a)

(b)

(c)

Figure 3.11: Surface electrostatic potentials of three proteins based on the molecular surface
after fractional PDE transform with α = 11.5 and t = 104 . (a) PDB code: 1VII; (b) PDB
code: 1FCA; (3) PDB code:1R69. Red: negative potential; blue: positive potential.

120

3.4.1.2

Electrostatic solvation free energy

With the solution of the Poisson equation, one can calculate electrostatic solvation free
energies of proteins by
∆Gp =

1
2

qj (φ(rj ) − φhomo (rj )),

(3.43)

j

where φ and φhomo correspond to the electrostatic potentials in the inhomogeneous and
homogeneous environments, respectively.
Electrostatic solvation free energies of twelve proteins are listed in Table 3.2. Results
are computed by using two types of surface deﬁnitions, namely, solvent excluded surface
generated by using the MSMS and the fractional PDE transform based surface. A good
consistence is obtained between these two methods.

3.4.2

Solvation analysis of 17 small compounds

Total solvation free energies, which include both polar (i.e., electrostatic solvation free energies ∆Gp ) and nonpolar solvation free energies Gnp , are measured in many experiments.
The latter component of solvation free energy can be modeled in many ways. One of the
most commonly used model can be expressed as [165, 37]

ρs U att dr,

Gnp = γArea + pVol +

(3.44)

Ωs

where γ is the surface tension, p is the hydrodynamic pressure, ρs is the solvent bulk density,
Ωs denotes the solvent accessible region, and U att (r) is the solute-solvent van der Waals
interaction potential at point r. The ﬁrst term γArea is the surface energy, which describes
the disruption of intermolecular and/or intramolecular bonds that occurs when the surface

121

of a molecule is created in the solvent. The second term pVol measures the mechanical work
of creating the vacuum of a biomolecular size in the solvent. The hydrophobic eﬀect in the
ﬁrst two terms are partially compensated by the third term Ωs ρs U att dr, which describes
the attractive dispersion eﬀects near the solvent-solute interface. The reader is referred to
Ref. [37] for detailed description of the van der Waals interaction.
One of the important applications of FPDE transform-based solute-solvent interface is
the calculation of total solvation energies. As stated in Section 3.3.4, all the terms in Eq.
(3.44) can be calculated from the molecular surfaces generated from FPDE transform of
various types of initial data, and according to Section 3.4.1.2, the analysis of electrostatic
free energies also requires the deﬁnition of solute-solvent interface.
Compound
glycerol triacetate
benzyl bromide
benzyl chloride
m-bis(triﬂuoromethyl)benzene
N,N-dimethyl-p-methoxybenzamide
N,N-4-trimethylbenzamide
bis-2-chloroethyl ether
1,1-diacetoxyethane
1,1-diethoxyethane
1,4-dioxane
diethyl propanedioate
dimethoxymethane
ethylene glycol diacetate
1,2-diethoxyethane
diethyl sulﬁde
phenyl formate
imidazole

Gnp ∆Gp
∆G
Exptl Error
3.38 -13.54 -10.16 -8.84 -1.32
2.00 -5.68 -3.67 -2.38 -1.29
1.97 -5.70 -3.73 -1.93
-1.8
4.38 -3.35
1.03
1.07
-0.04
2.85 -10.43 -7.57 -11.01 3.44
2.71 -8.56 -5.85 -9.76
3.91
2.09 -4.57 -2.48 -4.23
1.75
2.39 -8.77 -6.39 -4.97 -1.42
2.13 -4.52 -2.39 -3.28
0.89
1.38 -6.85 -5.47 -5.05 -0.42
2.67 -9.04 -6.37 -6.00 -0.37
1.34 -5.23 -3.89 -2.93 -0.96
2.31 -9.88 -7.57 -6.34 -1.23
2.14 -5.30 -3.15 -3.54
0.39
1.74 -2.93
-1.2
-1.43
0.23
1.98 -8.22 -6.25 -4.08 -2.17
1.11 -13.45 -12.34 -9.81 -2.53

Table 3.3: Computation of solvation free energies ∆G (kcal/mol) for the 17 small compounds
based on the molecular surface generated by fractional PDE transform with order α = 11.5
and propagation time t = 104 . Gnp : nonpolar component; ∆Gp : polar component; Exptl:
experimental data; Error: ∆Gp - Exptl.

122

In Table 3.3, there list calculated total solvation free energies of a set of 17 small compounds, compared with corresponding experimental data. The nonpolar free energy Gnp and
electrostatic energy ∆Gp are calculated by Eqs. (3.43) and (3.44), respectively, based on the
molecular surface based on the FPDE transform with order α = 11.5 and propagation time
t = 104 . The reader is referred to Ref. [38] for a detailed description of atomic coordinate and
charge parameters. The total solvation free energy ∆G = Gnp + ∆Gp . For comparison, the
corresponding experimental measurements (Exptl) of these compounds are also listed. We
can conclude that the application of FPDE transform-based molecular surface in calculations of free solvation energy is validated, because the agreement between the simulations and
experimental data is fairly good. The root mean square error (RMSE) of the computation
results is 1.77 kcal/mol and the average error is 1.42 kcal/mol. In comparison, Nicholls et al
[129] reported a RMSE of 1.87 kcal/mol for the same set of molecules by the linear Poisson
model and further, the RMSE was reduced to 1.71 ± 0.05 kcal/mol [129] by their explicit
solvent approach, which is much more expensive. Additionally, the RMSE in the current
work is only slightly greater than the one (RSME=1.76 kcal/mol) reported in [37], where the
biomolecular surface is deﬁned and generated from a nonlinear mean curvature ﬂow equation
so extra computational eﬀorts are necessary. Therefore, FPDE transform-based molecular
surface provides a relatively good foundation to predict the solvation energies for this set of
molecules.
Figure 3.12 depicts the surface morphologies and surface electrostatic potentials of 9
compounds. It is seen that fractional PDE transform oﬀers useful information of molecular morphology, which is crucial to the understanding of molecular physical, chemical and
biological properties.

123

3.5

Chapter conclusion remarks

The emergence of complexity in self-organizing biological systems poses exciting challenges
to their quantitative descriptions and predictions. Images, surfaces and visualization of
complex biomolecules, such as proteins, DNAs, RNAs, molecular motors and viruses, are
crucial to the understanding and conceptualization of biomolecular systems, which in turn
can have signiﬁcant impacts in biomedicine, rational drug design, drug discovery and gene
therapy. PDE transform is a new approach for mode decomposition of images, surfaces and
data [171, 169]. It makes use of arbitrarily high-order nonlinear PDEs for the control of
time-frequency localization and the regulation of spatial resolution. The PDE transform
has found its success in the analysis of non-stationary signals and noisy images [171, 169].
Recently, it has also been applied to the molecular generation of biomolecules [204]. However,
our previous PDE transform depends on the use of integer order PDEs of arbitrarily highorders. In this work, we extend the PDE transform to fractional PDE transform by using
fractional derivatives of arbitrarily high orders.
Fractional derivative or fractional calculus has received much attention in the past decade.
Its concept naturally arises in science and engineering, and includes integer derivatives as
special cases. However, most work in the ﬁeld involves only relatively low order fractional
derivatives. Fractional derivatives of orders higher than 2 are seldom used. It is not clear that
how to use arbitrarily high-order fractional derivatives in theoretical modeling and practical
computation.
To demonstrate the utility of fractional derivatives of arbitrarily high-orders, we propose high-order fractional PDEs and the corresponding fractional PDE transform. Using
the fractional variational principle, we construct nonlinear fractional PDEs based fractional

124

hyperdiﬀusion. Introducing an artiﬁcial time, the resulting high-order fractional PDEs are
converted to time-evolution fractional PDEs. Numerical techniques based on fast fractional Fourier transform (FFFT) is developed to compute the high-order fractional PDEs in
three-dimensional setting.
The proposed high-order fractional PDEs are applied to the surface construction of
macromolecular surfaces, which are crucial components in the implicit solvent models [67],
charge transport models [201] and variational multiscale models [180, 184]. We consider
three types of initial values to study the proposed high-order fractional PDEs. Additionally,
we examine the eﬀect of the orders of fractional PDEs on the surface morphology. It is found
that high-order fractional PDEs are crucial to the quality of molecular surfaces. We also test
the impact of the PDE integration time on surface generation. Moreover, we examine the
computational eﬃciency of the present method. Eﬃciency is one of major motivations for
developing new surface generation methods. It is found that the proposed high-order fractional PDEs are of linear scaling with respect to number of atoms in a molecule. We further
validate the present method by quantitative analysis of surface areas and surface enclosed
volumes of proteins. Finally, the surfaces constructed by the present approach is applied to
a couple of biophysical problems, namely, the electrostatic analysis via the Poisson equation
and the solvation analysis via a full solvation model. The results from these biophysical
problems indicate that the proposed high-order fractional PDEs are a robust and eﬃcient
method for macromolecular surface generation.

125

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3.12: Illustration of surface electrostatic potentials of nine small compounds at their
corresponding isosurfaces u(r, 104 ) = 0.5. (a) benzyl bromide; (b) bis-2-chloroethyl ether; (c)
1,1-diacetoxyethane; (d) diethyl propanedioate; (e) diethyl sulﬁde; (f) dimethoxymethane;
(g) m-bis benzene; (h) N,N-4-trimethylbenzamide; (i) phenyl formate.

126

Chapter 4
Modeling of Impact of geometric,
thermal and tunneling eﬀects on
nano-transistors

4.1

Introduction

Design, manufacture and operation of electronic devices, such as metal oxide semiconductor
ﬁeld eﬀect transistors (MOSFETs), have been down-scaled to nanometers, due to continuous
demand in improving the performance. The ultimate channel length of nano-MOSFETs is
expected to be below 10nm [65]. To this end, down-scaling of gate oxide, connecting material, doping size, operating voltages and other physical parameters is required [54, 43, 23].
At such a channel length, characteristics of devices quickly approaches the atomic scale and
the associated physical limits. Consequently, quantum mechanical eﬀects, including quantum coherence, correlation and/or entanglement, become signiﬁcant and will dramatically
impact macroscopic quantities of nano-devices, such as current-voltage characteristics and
conductance due to channel tunneling and gate leakage. These eﬀects may devastate classical
functions of MOSFETs. Recently, many properties of down-scaling associated devices are
investigated using quantum mechanical means. Some of these properties include (i) electron

127

scattering with semiconductor-on-insulator (SOI) interface roughness; (ii) electron scattering
under interactions with phonons; (iii) self-heating of devices; (iv) statistical ﬂuctuations due
to atomistic dopants; (v) choice of high-k materials for optimal device design and performances. Properties of nano-devices have been extensively studied for MOSFET, FinFET,
and various other silicon on insulator (SOI) devices in the literature [3, 84, 104, 189, 192, 193],
as well as two dimensional devices based on graphene [161].
Many tools in quantum mechanical theory have been developed in past decades in modeling and simulation of nano-scale devices. First of all, electronic structure in terms of
wavefunctions is required in most transport evaluations. However, quantum mechanical ﬁrst
principle description of nano-scale devices is intractable and approximations are indispensable. The coupled system including a Schr¨dinger equation and a Poisson equation usually
o
serves as one of the lowest levels of approaches [118, 160, 17, 100, 6, 59, 144, 12]. The
Schr¨dinger equation describes a single electron dynamics in a band structure of the solid,
o
whereas the Poisson equation provides the electrical potential of the device. Furthermore,
linear combination of bulk band (LCBB) method [97] , Hartree-Fock method [149] and the
density functional theory (DFT) [91, 107, 133] are theoretical models for many-body electronic interactions at diﬀerent levels of approximation accuracies and computational costs.
The non-equilibrium Green’s function (NEGF) method [168, 89, 45, 44] is currently the
main workhorse for electron structure and transport modeling of nano-devices. It is a general formalism to utilize the Fermi-Dirac statistics, and to include various interactions of
charge carries, such as scattering processes of particles (i.e., electrons and phonons), or
correlations due to the surroundings. On the other hand, some classical methods, such
as drift-diﬀusion models, Vlasov-Poisson equation, or Monte Carlo methods, can be modiﬁed by “quantum corrections” to study electron transport in the quantum ballistic regime
128

[88, 26, 48, 26, 4, 103]. Additionally, quantum Boltzmann equation, or known as WaldmannSnider equation [166, 150], can also give derivations of transport equations with quantum
corrections. The reader is referred to the literature [40, 60, 124, 84, 96, 35, 42] for many
other models and computational methods.
Recently, we have introduced a multiscale variational framework for nano-device modeling
[35]. In this framework, the macroscopic description of electrostatic potentials is integrated
with the microscopic description of electronic structure. We set up an energy functional for
the total energy of the system, which includes the electrostatic energy of the device and
the energy of the electron. The latter constitutes kinetic and potential energies. From the
variation of the energy functional with respect to electron wavefunction and electrostatic
potential, we derive coupled Kohn-Sham equation and Poisson equation for nano electronic
devices. When down to nano-scale, the semiconductor-insulator material interface can not
be considered to be smooth. We introduced interfaces to model the material separation
among diﬀerent components of device, which has an important impact on the governing
equations. For example, in order for the Poisson equation to describe the electrostatic
potential involving material interfaces, interface jump conditions are required to ensure the
well-posedness of the equation. Additionally, we introduce Dirac delta functions to model
random dopings. The resulting problems are elliptic equations with discontinuous coeﬃcients
and singular sources. A set of advanced numerical algorithms have been developed to improve
computational eﬃciency and accuracy in nano-device simulations: the matched interface and
boundary (MIB) method [207] and Dirichlet-to-Neumann Mapping (DNM) scheme [196]
ensure the accuracy of the electric ﬁeld with abruptly discontinuous dielectric permittivity
and individual dopants. We also used a mode-space decomposition method based on the
MOSFET geometry to reduce the fully 3D problem into the coupled 2D and 1D ones. A
129

self-consistent iteration method is used to solve the coupled equations.
The objective of the present work is to explore the impact of nano device design parameters on device the performance, including quantum tunneling eﬀects. We show that diﬀerent
conﬁgurations of channel geometries have impacts on device performance in terms of currentvoltage curves, or I-V curve, due to the geometry induced changes in electron conﬁnement
and electron-interface scattering. Additionally, only electron-potential ﬁeld scattering is
considered in the previous work, while in the formulation of current paper, electron-phonon
scattering is accounted as well. Further, temperature dependence and quantum tunneling
eﬀects in the channel length are tested for the whole system. In this study, we explore the
quantum tunneling eﬀect and macroscopic performances of a four-gate MOSFET model,
under diﬀerent channel geometry, phonon-electron interactions and temperature.
The rest of this paper is organized as follows. In Section 4.2, the theoretical model is
brieﬂy reviewed. A total free energy functional is set up and governing equations are derived
by the variational principle. Numerical algorithms for each governing equations and the
iteration procedure are described in Section 4.3. In Section 4.4, we design a number of interesting channel-crossing section shapes to understand the impact of geometry on a four-gate
MOSFET model with silicon/silicon-dioxide (Si/SiO2 ) interfaces. Numerical simulations are
also carried out to explore quantum tunneling eﬀects and phonon-electron interactions. The
paper ends with a conclusion in Section 4.5.
Our model is formulated with a uniﬁed free energy functional [35], which consists of
the free energy of electrons and electrostatic potential energy of the system. The whole
computational domain of the nano-device, denoted by Ω, is partitioned into two subdomains
due to the composited semiconductor and insulator materials. The silicon region and the
silicon dioxide part are denoted by ΩSi and ΩSiO2 , respectively. The semiconductor-insulator
130

interface separating these two subdomains is ΓSi/SiO , i.e., Ω = ΩSi ∪ ΩSiO2 and ΓSi/SiO =
2
2
ΩSi ∩ ΩSiO2 .

4.2

Theoretical models

Our model is formulated with a uniﬁed free energy functional [35], which consists of the free
energy of electrons and electrostatic potential energy of the system. The whole computational
domain of the nano-device, denoted by Ω, is partitioned into two subdomains due to the
composited semiconductor and insulator materials. The silicon region and the silicon dioxide
part are denoted by ΩSi and ΩSiO2 , respectively. The semiconductor-insulator interface
separating these two subdomains is ΓSi/SiO , i.e., Ω = ΩSi ∪ ΩSiO2 and ΓSi/SiO = ΩSi ∩
2
2
ΩSiO2 .

4.2.1

Free energy of electrons

For simplicity, we consider electrons as the major charge carriers in the device, i.e., density
and transport of holes are neglected. Due to the nano-scale of the device, electrons are
modeled quantum mechanically. The total Hamiltonian of the electrons includes kinetic
energy and potential energy, as described below.

4.2.1.1

Kinetic energy

Let Ψj be the Kohn-Sham orbitals of the electrons with energy Ej , then the kinetic energy
density of the electrons is
2 f (E
j
j

− µ)
| Ψj (r)|2 dr,
2m(r)

131

(4.1)

where m(r) is the eﬀective mass of electrons,

is the Planck constant, and µ is the Fermi

energy related to applied voltages. The general summation

includes the situation when
j

j is continuous, in which we have the identity in Rs for a function g(j)

g(j) =
j

2
(2π)2

g(j)(dj)s .

(4.2)

In Eq. (4.1), all energy levels Ej are weighted by the Fermi-Dirac function
1

f (Ej − µ) =
1 + exp

Ej −µ
kB T

,

(4.3)

where kB is the Boltzmann constant and T is the temperature.

4.2.1.2

Potential energy

The potential energy of electrons is associated with the density of electrons, which can be
represented by
|Ψj (r)|2 f (Ej − µ).

n(r) =

(4.4)

j

The potential is formulated in terms of the interactions of electrons among themselves and
with surroundings. First, electrons interact repulsively among themselves, then the energy
is
Uelec/elec =
where

1
2

q 2 n(r)n(r )
dr ,
|r − r |

(4.5)

is the dielectric permittivity.

Additionally, electrons interact with nuclei and the energy resulting from this interaction

132

is described as
Na

Uelec/nuclei = −
i

q 2 n(r)Zi
,
|r − ri |

(4.6)

where Zi q and ri (i = 1, 2, ..., Na ) are the charges and positions of the i-th nucleus, respectively.
Furthermore, devices are usually heavily doped, if we consider the speciﬁc N-type doping,
the interaction of electrons and dopants is
q 2 n(r)nd (r )
dr ,
|r − r |

Uelec/dop = −

(4.7)

where nd (r ) is the prescribed doping density. Note that at nano-scale, atomic dopants need
to be included. The formulation of interactions between electrons and individual dopants
can be accounted by the form of Eq. (4.6).
Since electrons are conﬁned in the channel with the semiconductor-insulator interface
ΓSi/SiO , we can write the total energy for electrons
2




Gelectron [n] =
Γ
Si/SiO2


+

1
2

2 f (E
j

− µ)
| Ψj (r)|2 + EXC [n] + Ephon [n]
2m(r)

q 2 n(r)n(r

j

)
dr −
(r)|r − r |

Na
i

q 2 n(r)Zi
(r)|r − ri |

−

q 2 n(r)n

d (r




)
dr dr,(4.8)

(r)|r − r |

where the exchange correlation EXC [n] and electron-phonon interaction Ephon [n] are also
included.

133

4.2.2

Electrostatic free energy of the system

The charge source of the nano-device system includes electrons, nuclei, and dopants, so
letting electrostatic potential of the system be u(r), we have the total free energy of electrics
as

Gelec
Γ

Si/SiO2

[u] =




Na

 (r)
2 − u(r) n(r)(−q) +
 dr,
| u(r)|
Zi qδ(r − ri ) + nd (r)q

 2
i

(4.9)
where (r) is the position-dependent dielectric permittivity since it is diﬀerent in the semiconNa

Zi qδ(r − r ).

ductor and the insulator regions. The charge density of nuclei is expressed by
i

Notice that u(r) can be expressed in a free space as

u(r) =

n(r )(−q)
dr +
(r)|r − r |

nd (r )q
dr +
(r)|r − r |

Na
i

Zi q
,
(r)|r − ri |

(4.10)

which is the formulation exactly representing the interactions between electrons and nuclei,
doping, and themselves in Eq. (4.8).

134

4.2.3

Total free energy functional

Combining Eqs. (4.8) and (4.9), and with relation (4.10), we obtain the total free energy
functional

GΓ

Si/SiO2

[u, n] =





2 f (E

− µ)
| Ψj (r)|2 + Uelec [n] + EXC [n] + Ephon [n]
2m(r)

j

j

Na

Zi δ(r − ri ) −

+ u(r)nd (r)q + u(r)q
i


+ λ n0 −
j

(r)
| u|2
2



f (Ej − µ)|Ψj (r)|2  dr.


(4.11)

The last line is the constraint for wavefunction Ψj (r), where λ is the Lagrange multiplier
and
1
n0 = √
2

mkB T 3/2
π 2

is the intrinsic density-of-state of an electron system.

4.2.4

The Poisson equation

Taking variation of Eq. (4.11) with respect to function u(r), we have the Poisson equation
δGΓ

Si/SiO2

δu

[u, n]
=0⇒−

· ( (r) u(r)) = qntotal ,

(4.12)

where
Na

ntotal = −n(r) + nd (r) +

Zi δ(r − ri ).

(4.13)

i

Since the permittivity (r) takes diﬀerent values in Si and SiO2 regions, it is discontinuous across the semiconductor-insulator interface ΓSi/SiO . We consider the following jump
2
135

conditions along ΓSi/SiO as
2

u+ (r) − u− (r) = 0,

[u] =
[

+

u · n] =

at ΓSi/SiO
2

u+ (r) · n − − u− (r) · n = 0,

(4.14)

at ΓSi/SiO ,

(4.15)

2

where n is the unit outer normal direction of the interface ΓSi/SiO , and the superscripts
2

“+” and “−” represent limiting values of the corresponding functions on ΩSi and ΩSiO2 ,
respectively. These conditions make the Poisson equation well-posed [196].

4.2.5

Generalized Kohn-Sham equation

Taking variation of Eq. (4.11) with respect to Ψ† (r), the complex conjugate of Ψ(r), we
have the generalized Kohn-Sham equation:

δGΓ [u, n]
=0⇒
δΨ†

2

−

·

2m

+ u(r)(−q) + UXC [n(r)] + Uphon [n(r)] Ψj (r) = Ej Ψj (r),
(4.16)

where UXC [n(r)] =

δEXC [n(r)]
δn

and Uphon =

δEphon [n(r)]
δn

are the exchange-correlation poten-

tial and electron-phonon potential, respectively.
The Poisson equation (4.12) and the Kohn-Sham equation (4.16) are to be solved iteratively to obtain the electrostatic potential and the wavefunction. This aspect is described in
the next section.

136

4.3

Computational algorithms

In the present paper, we neglect UXC [n(r)] for simplicity and leave it for future studies. In
the following discussion of numerical methods, we ﬁrst consider the model without phonon
eﬀect Uphon as in [35] and incorporate the phonon-electron interactions later on. As a result,
the equation (4.16) turns into the typical Schr¨dinger equation.
o
2

−

·

2m

+ u(r)(−q) Ψj (r) = Ej Ψj (r).

(4.17)

Note that Eq. (4.17) is coupled back to the Poisson equation (4.12), so it is nonlinear.

4.3.1

The Poisson equation

Individual dopants and their induced ﬂuctuations in device performance have been studied
elsewhere [35]. In the present work, we focus on the geometry eﬀects, tunneling behavior and
phonon-electron interactions, so only the continuous doping is considered and the Poisson
equation is taken as

−

·(

u) = −q(n + ND − NA ),

(4.18)

where ND and NA are given continuous doping density functions; the letter D and A
denotes donor and acceptor, respectively. For the boundary condition of Eq. (4.18), it takes
the Dirichlet boundary condition where external voltages are applied, denoted by Γvol , and
the Neumann boundary condition for the rest of the device.
Since ND and NA are given ﬁxed functions, it is convenient to decompose the potential

137

u as u = uﬁx + uloop . The potential component uﬁx solves the equation

(4.19)
(4.20)

uﬁx · n] =0 at ΓSi/SiO ,
2

·(

uﬁx ) = − q (ND − NA )
[uﬁx ] =0 at ΓSi/SiO

−

(4.21)

2

[

where the bracket [·] indicates the jump as deﬁned in Eqs. (4.14) and (4.15). We set the
boundary conditions

uﬁx =0 at Γvol

(4.22)

uﬁx · n =0 at other boundaries.

(4.23)

And the function uloop solves

−

·(

uloop ) = − qn(r)

(4.24)

[uloop ] =0 in ΓSi/SiO

(4.25)

uloop · n] =0 in ΓSi/SiO ,
2

(4.26)

2

[

with boundary conditions

uloop =uvoltage

at Γvol

uloop · n =0 at other boundaries.

(4.27)
(4.28)

The elliptic interface problems (4.19)-(4.23) and (4.24)-(4.28) are solved by the matched

138

interface and boundary (MIB) method [33, 207, 194, 196], which is a high-order numerical
scheme that can handle various challenges from material interface and have been successfully
applied to many other areas. For more details about the formulation of high-order MIB
methods and their applications, the reader is referred to Refs. [33, 207, 194, 196].

4.3.2

Scattering problem

Fully solving Eq.(4.16) in 3D is always a challenge in MOSFET simulations because of the
heavy computational cost. However, the computational complexity can be reduced by taking
advantage of the geometry and working principles of the devices. Particularly, if we consider
the direction connecting source and drain ends of the devices as the transport direction (denoted by x-direction), then electrons are conﬁned within y-z direction due to the insulator.
Therefore, the y-z directions are regarded as the conﬁned directions. Based on this considj

eration, the whole wave function can be written as Ψj (r) = Ψj,kx (r) = ϕj (y, z; x)ϕk (x),
x

and Eq. (4.16) is thus decomposed as, in the conﬁned y-z directions at each position x0 :
2

−

2my,z

d2
d2
+ 2
dy 2 dz

j

− qu(y, z; x0 ) ϕj (y, z; x0 ) = Escat (x0 )ϕj (y, z; x0 )

(4.29)

and along the transport x-direction,
d2
j
j
j j
+ Escat (x) ϕk (x) = Ek ϕk (x),
x
x x
2mx dx2
2

−

(4.30)

where mx and my,z are eﬀective electrons mass in the x and y-z directions, respectively.
When the electric ﬁled u(r) is available, we ﬁrst need to solve eigenvalue problems as Eq.
(4.29) at each position x0 , since the electrons are assumed not to penetrate the Si/SiO2 in-

139

j

terface. The resulting eigenvalues (or conﬁned energies) Escat (x0 ), j = 1, 2, 3, ... are discrete
and only the lowest several of them make signiﬁcant contributions due to the Fermi-Dirac
j

statistics, thus only the ﬁrst M functions Escat (x), j = 1, 2, ..., M and corresponding eigenj

functions ϕk (x) are typically utilized. Here M is usually less than 5 due to the fast increase
x

j

in energy Escat (x) as j increases.
j

Equation (4.30) is a scattering problem so that the eigenvalue Ek is continuous; the
x

function

j
Escat (x)

obtained from Eq. (4.29) is sometimes called subband energy and considj

ered as transport potential energy. The wavefunction ϕk (x) is not directly solved, in stead,
x

the Hamiltonian of electrons in scattering direction, Hscat , is given as

Hscat = −

d2
j
+ Escat (x),
2mx dx2
2

(4.31)

and the j-th scattering electron density in Eq. (4.33) is calculated as [44, 35]:

j

nscat (x) =

1
j
j
j
f (Hscat − µS I)AS (Ek ) + f (Hscat − µD I)AD (Ek ) dEk ,
x
x
x
2π R

(4.32)

j

where I is the identity operator, AS (Ek ) and µS are the spectral operator and Fermi-level
x

j

at the source end, respectively. Similar notations, AD (Ek ) and µD are used for the drain
x

end.
By this means, the total electron density is given by
M
j

j

|ψj (y, z; x)|2 |ψk (x)|2 f (Ek − µ)

n(r) =

x

x

j
M
j

|ψj (y, z; x)|2 nscat (x),

=
j

140

(4.33)

and total electron current is computed in the formula of
∞

Itot = −

q

−∞

j

j

TrTDS (Ek )[f (Htot − µS I) + f (Htot − µD I)]dEk ,
x

x

(4.34)

j

where Htot = Hscat + Ek is the total Hamiltonian, Tr(·) is the trace operator, and TDS is
x

the transmission operator [44, 35].

4.3.3

Modeling the phonon-electron interaction

Acoustic and optical phonon scattering is an important supplementary component of the
interactions of electron with the whole environment. In the Boltzmann transport equationbased classical models [41], the phonon-electron interaction appears in the collision operators
of the Boltzmann equation, while in the non-equilibrium Green’s function (NEGF) model
for nano-devices [44], this interaction is considered as a self-energy function. Equation (4.11)
has taken into account the phonon-electron interaction in the total energy framework and
our goal here is to develop a direct formulation.
First we assume that, as in [41], all the phonons have a single frequency ωp and the
number density of phonons is homogeneous in the device. Thus we can consider the phononelectron interaction mainly in the transport direction, then the Eq. (4.30) is modiﬁed by
incorporating Uphon :
d2
j
j
j j
+ Escat (x) + Uphon ϕk (x) = Ex ϕk (x).
x
x
2mx dx2
2

−

(4.35)

Further, it is assumed that the phonon density depends on its energy Ep = ωp and

141

follows the Bose-Einstein statistics,

Np = exp

Ep
KB T

−1

−1

.

(4.36)

Therefore, we have Uphon
Uphon = Np Ep .

(4.37)

In the present model, two statistics, the Bose-Einstein statistics for phonons and the FermiDirac statistics for electrons, are involved. Since both statistics depend on temperature, T
will play an important role in the model performance.

4.3.4

Self-consistent iteration

The self-consistence of Eqs. (4.12) and (4.16) are numerically achieved by the Gummel-like
iteration [47, 35] and the procedure is outlined as following:
• Step 0 (Solution of uﬁx ) Solving Eq. (4.19-4.23) is out of the whole iteration loop and
once uﬁx is obtained, it remains ﬁxed in every step in the proceeding iterations;
loop

• Step 1 (Solution of uloop by inner iteration): We assume that nk and uk

,k =

0, 1, 2, ... are available from the previous steps, then solve the following non-linear
loop

Poisson equation for uk+1 :

−

·(

loop

uk+1 ) = −qn0 F1/2 

142

(uﬁx

loop
+ uk+1

kB T

− ζk )q


,

(4.38)

where ζk is the quasi-Fermi potential in the k-th step calculated through

loop

ζk (r) = uﬁx + uk

−1
+ F1/2 (nk /n0 ),

(4.39)

and F1/2 is the Fermi-Dirac integral of order 1/2 deﬁned as
∞

F1/2 (x) =

0

y 1/2 dy
.
1 + ey−x

(4.40)

The non-linear Poisson equation (4.38) is iteratively solved with the corresponding
interface and boundary conditions in Eqs. (4.25)-(4.28).
• Step 2 (Eigenvalue problem): We solve Eq. (4.29) in the y-z directions for ψj (y, z; xi )
j

and Escat (xi ) along the x-direction at each grid point xi .
• Step 3 (Scattering problem and electron density): We solve the scattering problem
(4.30)-(4.32); using (4.33), we obtain the electron density nk+1 (r) in the k + 1-th step.
• Step 4 (Convergence check): The convergence is checked by ||nk+1 − nk || < δ, where
δ is a given criteria. If convergence is achieved, then go to Step 5, otherwise go back
to Step 1.
• Step 5 (Electronic current): After we obtained the convergent potential and electron
density, we can obtain the electronic current by Eq. (4.34).
Figure 4.1 gives an illustration of the workﬂow for the current self-consistent iteration scheme.

143

Solve for uﬁx
Solve uloop by inner iteration

Solve eigenvalue problem

Solve scattering problem

||nk+1 − nk || < δ ?

Outer iteration

no

yes

Output electronic current

Figure 4.1: The ﬂowchart of Gummel-like self-consistent iteration.

4.4
4.4.1

Numerical experiments and discussions
A four-gate MOSFET model

To study the channel geometric eﬀects, electron tunneling and electron-phonon interactions
by the proposed model, we consider a simple four-gate MOSFET device, i.e., silicon nanowire
transistor and present our numerical results for the proposed models. Figure. 4.2(a) illustrates the overall conﬁguration of this type of MOSFET: the semiconductor channel (silicon)
in cylinder shape is bounded by the insulator material (silicon dioxide). Functions of the
nano-devices are manipulated by two types of voltages: the source-drain denoted by VDS
and the voltage applied on the all-around gate denoted by VGate . Without loss of generality,
we set x-direction as the electron transport direction and y-z directions as the conﬁned di144

Figure 4.2: Illustration of a four-gate MOSFET with a cylinder channel: (a) Overall 3D
conﬁguration; (b) View of the x-z cross section.

rections. Figure 4.2(b) gives a detailed demonstration of the device in the x-z cross-section,
from which it can be seen that N-type doping (denoted by N) is applied to the two ends of
the channel. In the current work, we also apply the present model to channels with diﬀerent
shapes of y-z cross-sections other than the disk. As shown in Figure 4.3(a) and (b), the
channels with y-z cross-sections of a circle and a ﬂower-like geometry, are also considered.
In the sequel, the three types of channels with diﬀerent cross-sections in y-z directions are
called the cube, cylinder and ﬂower-like channels as explained in Figure 4.2. The total length
of the device along x-direction is 30nm with 10nm for each doping area and channel.
For computational convenience, the x-z cross-section of the MOSFET shown in Figure
4.2 (b) is simpliﬁed as the one in Figure 4.4. The channel (LGHK) is sandwiched by the

145

5

z (nm)

4
3
2
1
0
0

1

2
3
y (nm)

4

5

4

5

(a)

5

z (nm)

4
3
2
1
0
0

1

2
3
y (nm)
(b)

Figure 4.3: (a) y-z cross-section of a channel: circle; (b) y-z cross-section of channel: ﬂowerlike shape.

insulator layers ABGL and KHCD, and LG and KH are the two Si/SiO2 interfaces. The
boundary segments EF and JI are considered where the gate voltages are applied while LK
and GH are regarded as source and drain ends, respectively. Hence, Figure 4.4 also indicates
the boundary conditions of the Poisson equation: Dirichlet boundary conditions VGate are
applied on EF and JI. Similarly, Dirichlet boundary conditions VDS are applied on LK and

146

GH. Homogeneous Neumann boundary conditions are imposed on the rest of the boundary.
In solving the Schr¨dinger equation, electrons are assumed to be conﬁned by LG and
o
KH, along with the boundaries of diﬀerent shapes in the y-z cross-section (n(r) = 0 at these
boundaries), while scattering boundary conditions are used at LK and GH.

Figure 4.4: Demonstration of the simpliﬁed computational domain for the MOSFET. Dark:
insulator layer; grey: doped channel region; and white: un-doped channel region. Gate
voltages are applied at boundary segments EF and IJ, and source-drain voltages are applied
on the boundary segments GH and LK.

In the present model, the parameters for the equations are set as follows. For the Poisson
equation, the dielectric constants for silicon layer and silicon dioxide layer are Si = 11.7 0
and SiO2 = 3.9 0 , respectively, where 0 is the dielectric constant for the vacuum. The
reference continuous N-doping is taken as ND (r) = 2 × 1020 /cm3 , NA (r) = −1 × 1020 /cm3 .
For the eigenvalue problem and scattering problem, the ﬁrst three subbands are adopted since
they dominantly contribute to electron density and current over higher ones. The electron
eﬀective mass is taken as mx = 0.50m0 in the transport direction and myz = 0.318m0 in
the conﬁned directions. Room temperature of T = 300K is assumed except when speciﬁed.

147

4.4.2

General results of electric ﬁeld and electron density

−0.4

eV

−0.6
−0.8
−1
4
30

3

20

2
nm

1 0

10
nm

1e26 /m3

(a)

2
1
0
4
30

3

20

2
nm

1 0

10
nm

(b)
Figure 4.5: Electrostatic and electron density proﬁles from the x-z cross-section view at
y = 2.5nm in the silicon area of the MOSFET with the cylinder channel (Source-drain
voltage VDS = 0.4V and gate voltage VGate = 0.4V). (a) Electrostatic potential energy
(eV); (b) Electron density (1e26/m3 ).

148

We ﬁrst present some fundamental proﬁles of electric ﬁeld and electron density for the
MOSFET with a cylinder channel. In order to give a clear description, 2D data in the x-z
cross-section (at y = 2.5nm) of the silicon area of the MOSFET are shown in Figs. 4.5(a)
and (b), for electric ﬁeld and electron density, respectively. In Figure 4.5(a), the electrostatic
potential energy of electrons is relatively lower in the source and drain regions due to the
dopants while remains as an energy barrier in the middle. Electrons are injected from the
source and collected at the drain; they have to overcome the potential barrier in the middle
in order to transport along the devices. Correspondingly, as shown by Figure 4.5(b), the
resulting electron density is lower in the middle and higher in the source and drain regions.
Since we assume that electrons cannot penetrate the semi-conductor-insulator interface, the
electron density is zero in the silicon dioxide region, which is not displayed in the ﬁgures.
In there ﬁgures, the source-drain voltage diﬀerence is taken as VDS =0.4V and the gate
voltage is taken as VGate = 0.4V. Later on we will see the performance of the device can be
controlled by adjusting VDS and VGate .
The potential barrier in the channel can also be illustrated by the subband energy proﬁles.
Due to the nano scale of the devices, there are strong conﬁnement eﬀects in the y-z crosssection, thus the total transporting energy can be approximated as discrete subband energies.
Because of the Femi-Dirac distribution, only the ﬁrst or lowest several of them have signiﬁcant
contributions to the total electron current. Figure 4.6 illustrates the ﬁrst subband energy
proﬁles of the MOSFET at two diﬀerent gate voltages, with the source-drain voltages ﬁxed
at VDS = 0.4V. It can be concluded from the ﬁgure that increasing the gate voltage will
lower down the potential barrier in the device channel and thus lead to a larger electron
current, which can be further seen in the I-VGate curve in Figure 4.7 (b).
Figure 4.7 illustrates the impacts of gate voltage and source-drain voltage diﬀerences on
149

Subband Energy(eV)

0.5

0

−0.5
0

10
20
Transport Direction(nm)

30

Figure 4.6: Subband energy proﬁle along the transport direction under gate voltages VGate =
0.1V(Dashed line) and VGate = 0.4V(Solid line) for the MOSFET with the cylinder channel
(source-drain voltages is ﬁxed at VDS = 0.4V).

4

3

3

I (µ A)

5

4
I (µ A)

5

2
1
0
0

2
1

0.1

0.2
VDS (V)

0.3

0
0

0.4

(a)

0.1

0.2
VGate ( V)

0.3

0.4

(b)

Figure 4.7: Current-voltage (I - V) curves of the MOSFET with the cylinder channel. (a)
I-VDS curve (gate voltage ﬁxed at VGate =0.4V); (b) I-VGate curve (source-drain voltage
ﬁxed at VDS =0.4V).

the electron current, respectively. It is observed from Figure 4.7(a) that, when the gate
voltage is ﬁxed at VGate =0.4V, the electron current increases with the source-drain voltage

150

5

60
50

4
I (µA)

z (nm)

40
3
2

outside circle
square
inside circle
flower

30
20

1

10

0
0

1

2
3
y (nm)

4

0
0

5

(a)

0.1

0.2
VGate (V)

0.3

0.4

(b)

Figure 4.8: I-VGate curves of MOSFETs with diﬀerent geometries of channel (VDS =0.4V).
(a) Four types of cross-sections channels; (b) Corresponding I-VGate curves.

diﬀerence and achieves a saturation later on. In contrast, when the source-drain voltage is
ﬁxed at VDS =0.4V, the electron current grows dramatically when the gate voltage VGate
passes a certain threshold, as shown in Figure 4.7 (b).

4.4.3

Performances of nano-MOSFETs with diﬀerent channel geometries

Performances of nano-MOSFEs with diﬀerent channels (cube, cylinder, and ﬂower-like channels) are examined in terms of I-V curves.
In Figure 4.8 we compare the I-V curves of the channels that have diﬀerent areas of the
cross-sections. Four types of channel cross sections, including a small circle, a large circle,
a square, and a ﬂower-like shape, are tested in our experiments. These silicon channels are
embedded in the insulator materials (SiO2 ) of the same dimension (5nm by 5nm). The
interfaces enclose diﬀerent silicon areas and their relations are shown in Figure 4.8(a). In

151

4

3

3

I (µ A)

5

4
z (nm)

5

2

2
1

1
0
0

circle
square
flower

1

2
3
y (nm)

4

0
0

5

(a)

0.1

0.2
VGate (V)

0.3

0.4

(b)

Figure 4.9: The I-VGate curves of MOSFETs with diﬀerent channel cross-sections but with
the same cross-section area ( source-drain voltage ﬁxed at VDS =0.4V).

Figure 4.8(b), the I-V curves of the nano-MOSFETs with the corresponding y-z cross-sections
are displayed. It can be easily concluded that given the same gate and source/drain voltage,
the smaller of channel cross-section area is, the lower electron current is obtained in the
corresponding MOSFET.
It is also interesting to explore the performances of the nano-MOSFETs with diﬀerent
channel geometries but with the same area of the y-z channel cross-section. This set of
experiments are designed to test geometry eﬀects of the Si/SiO2 interfaces. We still consider
the nano-MOSFETs with square, circle and ﬂower-like channel cross sections and examine
their corresponding I-V curves in Figure 4.9. The areas enclosed by the square, circle and
ﬂower-like curve, see Figure 4.9 (a), are the same as the small circle does in Figure 4.8(a).
Comparing the three solid I-V curves in Figure 4.8(b) and Figure 4.9(b), respectively, we
ﬁnd that although the areas of the cross-sections are the same, the electron currents at
various gate voltages from the cylinder channel are the highest, then followed by the currents

152

from channels with square and ﬂower-like cross sections. Therefore, we conjecture that for
the cross-sections of channels in the conﬁned directions, the smaller perimeter of Si/SiO2
interface will result in larger magnitude of electron current, when the area are the same.
The performances of the three types of channels are further analyzed in terms of their
corresponding subband energies along the transport direction. In Figure 4.10, the subband
energy proﬁles of the channels with the same area of cross-sections but with square, circle
and ﬂower-like cross section shapes, are plotted along the conﬁned direction. As can be
seen in Figure 4.10 (a), the barriers of subband energies (i.e., eigenvalues) of the channels
with ﬂower-like and circle cross sections are the highest and lowest, respectively, which are
consistent with the corresponding the smallest and largest magnitudes of electron currents.
The conﬁnement eﬀects of the three diﬀerent cross-sections with the same area are shown in
Figure 4.10 (b), in which the diﬀerence between the ﬁrst (u1 (x)) and the second (u2 (x))
scat
scat
subband energies is plotted for each channel. The diﬀerences between the ﬁrst two subband
energies from the channels with ﬂower-like and circle cross sections are the largest and
smallest, respectively, which indicate that the conﬁnement of the channels with ﬂower-like
and circle Si/SO2 interfaces are the strongest and weakest among the three geometries,
although their enclosed areas are the same.

4.4.4

Quantum tunneling eﬀects

A signiﬁcant diﬀerence between nano-devices and classical devices is the possible electron
tunneling eﬀects. From Figure 4.5(a) or 4.6, we can see that there is an electric energy
barrier in the channel. We denote the maximum of subband energy by Emax , i.e., Emax =
j

sup Escat (x0 ). In classical models, any electron with energy lower than Emax from the source
x0

can not overcome the barrier and arrive the drain. However, in a quantum mechanical model,
153

Subband Energy(eV)

0.2

circle
square
flower

0
−0.2
−0.4
−0.6
−0.8
0

5

10
15
20
Transport Direction(nm)

25

30

Subband Energy differences(eV)

(a)
0.8

circle
square
flower

0.75
0.7
0.65
0.6
0.55
0.5
0

5

10
15
20
Transport Direction(nm)

25

30

(b)
Figure 4.10: Subband energy proﬁles along the transport direction for MOSFETs with different geometries of channel cross-sections(source-drain voltages is ﬁxed at VDS = 0.4V and
gate voltages VGate = 0.4V). (a) First subband energies of diﬀerent channels; (b) Diﬀerence
of ﬁrst two subband energies of each channel.

154

electrons with energy less than Emax may be able to pass the barrier and result in electron
current. This electron tunneling current cannot be neglected in the nano-scale devices.
A method to study the quantum tunneling is to investigate the contribution to the
current density by electrons of various energies. The theoretical total current density Itot
is computed by the formula in Eq. (4.34). For convenience, we only consider the transport
energy E and use the formulation:
∞

Itot = −

q

−∞

TrTDS (E)[f (Htot − µS I) + f (Htot − µD I)]dE,

(4.41)

deﬁne the total current contributed by the tunneling eﬀects, Iqt , by
Emax

Iqt = −

−∞

q

TrTDS (E)[f (Htot − µS I) + f (Htot − µD I)]dE,

(4.42)

and denote the ratio of quantum tunneling current by

Rqt =

Iqt
.
Itot

(4.43)

Obviously, if Iqt = 0, the mechanism is entirely classical. By contrast, in the present quantum
model, Iqt > 0 because of the quantum tunneling eﬀect, i.e., there is a nonzero probability for
the electron with energy lower than the potential barrier to penetrate through the channel.
In Figure 4.11(a), we demonstrate the electron current (I) as a function of electron energy
(E) for a nano-MOSFET with cylinder channel under VGate = 0V and VDS = 0.4V. It can
be concluded that the electron current increases with higher electron energy (this is quite
understandable because the electron with a higher energy can overcome the potential barrier
more easily), and saturates fast when the energy is over around E = 0.5V. This is due to the
155

Current (1e−7 µ A)

12
10
8
6
4
2
0
0

0.2

0.4

0.6

0.8

1

E
(a)

Normalized current

1
0.8
0.6
0.4
0.2
0
0

0.1

0.2

0.3

0.4

0.5

E
(b)
Figure 4.11: Relation between electron current (µA) and electron energy E (eV) in the MOSFET with the cylinder channel. (a)Gate voltage and source-drain voltage are respectively
VGate = 0.0V and VDS = 0.4V. The red dot indicates the current w.hen electron energies
equal to the maximum energy barrier Emax (horizontal coordinate). (b)under various gate
voltages VGate (source-drain voltages ﬁxed at VDS = 0.4V). The red dots from right to left
correspond to VGate = 0, 0.05, 0.1, ..., 0.4V.

156

VGate (V)
0
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40

tunneling current(µA)
1.92 × 10−7
1.36 × 10−6
9.57 × 10−6
6.61 × 10−5
4.44 × 10−4
2.73 × 10−3
1.45 × 10−2
5.70 × 10−2
2.26 × 10−1

total current(µA)
1.20 × 10−6
8.79 × 10−6
6.42 × 10−5
4.69 × 10−4
3.40 × 10−3
2.43 × 10−2
1.66 × 10−1
9.60 × 10−1
4.07

tunneling ratio
15.94%
15.48%
14.92%
14.11%
13.04%
11.25%
8.78%
5.94%
5.55%

Table 4.1: Tunneling ratio for the cylinder channel at diﬀerent VGate (VDS =0.4V).

Fermi-Dirac distribution; the probability of electrons with higher energies is lower so they
do not contribute signiﬁcantly. In this simulation, the maximum energy barrier Emax of the
ﬁrst subband is Emax = 0.388eV; the point which has this value as horizontal coordinate on
the I-E curve is marked with a red dot. The tunneling current Iqt in this circumstance is
apparent since the electron current I is non-zero (up to 1.92 × 10−2 µA) when E < Emax .
From the total current Itot = 1.20 × 10−6 µA, we compute Rqt =15.94%.
The tunneling current Iqt and ratio of tunneling current Rqt depend on the gate voltage
of the nano devices. In Figure 4.11(b) there collects the similar current-energy curves as in
Figure 4.11(a), but with VGate takes values from 0.0V to 0.4V. Since the currents are significantly diﬀerent under values of VGate , these curves are normalized for better comparison.
In this sense the vertical coordinate of the red dots exactly represents the tunneling ratio
Rqt . We see that as VGate varies from 0.0V to 0.4V, the I-E curves shift to the left. This is
consistent with the fact that the electron with the same energy are easier to overcome the
potential barrier with assistance of voltage gate. It is also observed that the tunneling ratio
decreases as VGate increases. The detailed tunneling currents and tunneling ratios under
diﬀerent VGate are listed in Table 4.1.
157

0.3

flower
circle

Tunnelling Ratio

0.25
0.2
0.15
0.1
0.05
0
0

0.1

0.2
VGate (V)

0.3

0.4

Figure 4.12: Relation of electrons tunneling ratio and gate voltages VGate for nanoMOSFETs with cylinder and ﬂower-like channels (source-drain voltage ﬁxed at VDS = 0.4V).

Figure 4.12 illustrates the tunneling ratios for the two diﬀerent channels, the cylinder
channel and the ﬂower-like channel, whose y-z cross-sections are shown in Figure 4.8 (a). For
both devices, the tunneling ratio decreases as gate voltages VGate increases. Furthermore,
the tunneling ratio of the device with ﬂower-like channel is higher than the one with cylinder
channel. This is the evidence that stronger conﬁnement of channels will result in more
signiﬁcant quantum eﬀects.

4.4.5

Phonon-electron interactions

In this section, we investigate the phonon eﬀects based on the modiﬁed scattering equation
in Eq. (4.35). As seen from Eqs. (4.35)-(4.36), the potential energy component in Eq. (4.35)
will be higher than that in Eq. (4.30). There are essentially a couple of factors to aﬀect the
Uphon , the phonon energy and phonon number. In the present work, the phonon energy Ep

158

is assumed to take a constant in each simulation.
Figure 4.13(a) gives a comparison of I-VGate between the models with and without
phonon-electron (p-e) interactions. Compared to the model without p-e interactions, the
model with p-e interactions gives smaller electron currents. Moreover, the impact in the
case of Ep = 0.03eV is greater than that of Ep = 0.063eV. The impact of Uphon can be
seen as a trade-oﬀ between phonon energy and phonon number. In case of higher phonon
energy, there are a smaller number of phonons based on Bose-Einstein statistics in Eq.
(4.36). Eventually, the magnitude of Uphon decreases as Ep increases from Ep = 0.03eV
to Ep = 0.063eV, which results in a smaller impact for scattering potential barrier in the
scattering problem and furthermore higher current. Figure 4.13(b) gives results for the tunneling ratio. Although the tunneling ratio keeps decreasing as I-VGate increases, there is no
signiﬁcant diﬀerence between models with and without p-e interactions.

I (µ A)

4

0.16

without p−e interaction
Ep=0.063eV
Ep=0.03eV

0.14
Tunnelling Ratio

5

3
2
1
0
0

0.12
0.1
0.08
0.06

0.1

0.2
VGate (V)

0.3

0.04
0

0.4

(a)

without p−e interaction
Ep=0.063eV
Ep=0.03eV
0.1

0.2
VGate (V)

0.3

0.4

(b)

Figure 4.13: Eﬀects of phonon-electron interactions in I-V curves and electron tunneling ratio
for the cylinder channel. (a) I-VGate curves with diﬀerent electron-phonon interactions; (b)
Relations of electron tunneling ratios and gate voltages VGate with diﬀerent electron-phonon
interactions (source-drain volate ﬁxed at VDS =0.4V).

159

Figure 4.14 illustrates the impact of phonon eﬀects at diﬀerent temperatures. The result
in Figure 4.14(a) shows that current density increases with the increase of temperature. It
is consistent with the fact that the electronic conductivity of semiconductor enhances as the
temperature increases.
Figure 4.13(b) suggests that the tunneling ratio decreases as the temperature increases. This
trend is consistent with the quantum mechanical principle that quantum eﬀects play a more
important role under lower temperature.
3

0.6

without p−e interaction
Ep=0.063eV
Ep=0.03eV

0.5
Tunnelling Ratio

2.5
2
I (µ A)

without p−e interaction
Ep=0.063eV
Ep=0.03eV

1.5
1
0.5

0.4
0.3
0.2
0.1

0
100

200

300
400
T (K)

500

0
100

600

(a)

200

300
400
T (K)

500

600

(b)

Figure 4.14: Temperature eﬀects on electron currents under diﬀerent phonon-electron interactions for the cylinder channel (VDS =0.4V, VGate =0.3V). (a) Current-temperature curves;
(b) Electron tunneling ratios under diﬀerent temperatures.

4.5

Chapter conclusion remarks

We present mathematical modeling and numerical simulations of semiconductor-insulator
geometric eﬀects, quantum tunneling currents and phonon-electron interactions for nanoscale electronic devices, using a four-gate MOSFET model. Transistors in nano-scales will

160

not function in analog with conventional devices by classical physics and will lead to new
breakthroughs in devices science and technologies. For example, channel tunneling or gate
leakage under various voltage settings of ultimate MOSFETs and other nano-transistors
will be of main concerns of the modeling and simulations. Other factors may aﬀect the
performances of nano-MOSFETs at such a scale or under quantum mechanics could include
channel surface roughness, and phonon scattering eﬀects. However, quantum simulations are
intractable for individual particles because of the enormous degrees of freedom. We proposed
a two-scale energy functional which describes the electrons and the continuum electrostatic
potential of the nano-electronic device [35]. In our framework, the microscopic description
of electrons and the macroscopic description of the continuum electrostatic potential are
integrated on an equal footing at nano-scale.
The present work, based on our framework of free energy functional, addresses the modeling issues of various types of Si/SiO2 geometries, tunneling ratios and phonon-electron
interface eﬀects. Material interfaces are integrated in the energy formulation and serve as
continuity conditions of the Poisson equation and boundary conditions of the Schr¨dinger
o
equation. Thus, electric ﬁeld, physical conﬁnement of electrons, and in further macroscopic
performances of nano-MOSFETs will depend on these microscopic details. Phonon-electron
interactions are modeled as part of transport energy in a fashion of density functional theory,
in which density of phonons are assumed to be homogeneous in the channel and follows the
Bose-Einstein statistics. Finally, the tunneling eﬀects of devices are measured as tunneling
ratios of electron currents.
A four-gate MOSFET with diﬀerent shapes of channels is considered in the three-dimensional
numerical simulation in the present work. In our computations, advanced interface-problem
solvers are utilized to solve the Poisson equation for the electrostatic potentials with three
161

semiconductor-insulator geometries. The Schr¨dinger equation is solved based on the physo
ical properties that electrons have continuous and discrete spectrum along transport and
conﬁned direction, respectively. The self-consistence of the coupled system is achieved by
the so-called Gummel-like iteration. The performances of the nano-device are examined
in terms of current-voltage curves with three types of Si/SiO2 interfaces, a wide range of
source-drain and gate voltages, and diﬀerent phonon energies.

162

Chapter 5
Thesis achievements and future work

5.1

Contributions

The main contributions of this thesis are as follows. First, we propose a new electrostatic solvation free energy functional to take into consideration the eﬀect of hyperpolarizations. The
nonlinear Poisson equation(NLPE) is derived from electrostatic solvation free energy functional based on the variational principle. Compared to the classic Poisson equation, which
is derived based on electric polarizations in a linear, isotropic, and homogeneous dielectric
medium, our NLPE partially accounts for the hyperpolarizations, which are important under a strong electrical ﬁeld, or involving highly charged nonlinear inhomogeneous media.
Moreover, our nonlinear Poisson equation gives rise to a smooth dielectric function, which
overcomes computational diﬃculty of employing a sharp solvent-solute boundary. A nonpolar solvation model is further formulated on the basis of the dielectric proﬁle obtained from
the nonlinear Poisson model by using geometric measure theory. The proposed solvation
models are extensively validated with the Kirkwood model and experimental measurements
of 17 molecules. Our new solvation models out-perform the classic Poisson equation based
solvation model. In fact, its performance is very close to the explicit solvation model, which is
much more computationally expensive. Application of the proposed nonlinear Poisson model
is considered to electrostatic analysis of 20 proteins. In addition, the test results for the set
of 21 compounds at diﬀerent temperatures further validate that the proposed NLPE model
163

does a good job in modeling the temperature eﬀect on the solvation free energy. The good
agreement between our results and experimental data as well as theoretical results suggests
that the proposed nonlinear Poisson model is a potentially powerful model for electrostatic
analysis involving hyperpolarization eﬀects.
Second, we employ the PDE transform to eﬃciently suppress Gibbs’ oscillations in the
numerical solution of hyperbolic conservation laws. Speciﬁcally, during the time integration
of a hyperbolic conservation law system, an intermediate numerical solution at a given time
step is used as an initial data for a special evolution PDE. Then the solution of such an
evolution PDE is accepted as an updated numerical solution at the given time step. Two
techniques are proposed to improve the eﬃciency of the present approach. First, we use an
adaptive measure of total variations to automatically determine whether the PDE transform
is needed at each time step. Additionally, we utilize a fast PDE transform, which oﬀers the
analytical solution of an arbitrarily high order evolution PDE in the Fourier representation.
This technique bypasses the stability constraint of solving high order evolution PDEs. We
have employed a variety of benchmark in the present work to validate the proposed approach,
ranging from scalar conservation law systems to Euler equations in one and two spatial
dimensions. Among these problems, some typically prefer low-order shock capturing schemes;
whereas others are well known to require high order numerical methods. For example, low
order schemes will severely damp the amplitude of the entropy waves in the shock-entropy
interaction described by the Euler equation. The proposed PDE transform provides some
of the best results for solving the Burgers’ equation with non-convex ﬂux. Furthermore,
only about 5 points per wavelength (PPW) is needed for the present approach to handle the
interaction of shock-entropy waves and shock-vortex interactions, which is among the fewest
number of PPW to out knowledge so far.
164

Third, we propose high-order fractional PDEs and the corresponding fractional PDE
transform. Using the fractional variational principle, we construct nonlinear fractional PDEs
based on fractional hyperdiﬀusion. Introducing an artiﬁcial time, the resulting high-order
fractional PDEs are converted to time-evolution fractional PDEs. Numerical techniques
based on fast fractional Fourier transform (FFFT) is developed to compute the high-order
fractional PDEs in three-dimensional setting. The proposed high-order fractional PDEs are
applied to the surface construction of macromolecular surfaces, which are crucial components in the implicit solvent models. We consider three types of initial values to study the
proposed high-order fractional PDEs. Additionally, we examine the eﬀect of the orders of
fractional PDEs on the surface morphology. It is found that high-order fractional PDEs are
crucial to the quality of molecular surfaces. We also test the impact of the PDE integration time on surface generation. Moreover, we examine the computational eﬃciency of the
present method. Eﬃciency is one of major motivations for developing new surface generation methods. It is found that the proposed high-order fractional PDEs are of linear scaling
with respect to number of atoms in a molecule. We further validate the present method by
quantitative analysis of surface areas and surface enclosed volumes of proteins. Finally, the
surfaces constructed by the present approach is applied to a couple of biophysical problems,
namely, the electrostatic analysis via the Poisson equation and the solvation analysis via a
full solvation model. The results from these biophysical problems indicate that the proposed
high-order fractional PDEs are a robust and eﬃcient method for macromolecular surface
generation.
Finally, we present mathematical modeling and numerical simulations of semiconductorinsulator geometric eﬀects, quantum tunneling currents and phonon-electron interactions
for nano-scale electronic devices, using a four-gate MOSFET model. The present work,
165

based on our framework of free energy functional, addresses the modeling issues of various
types of Si/SiO2 geometries, tunneling ratios and phonon-electron interface eﬀects. Material interfaces are integrated in the energy formulation and serve as continuity conditions
of the Poisson equation and boundary conditions of the Schr¨dinger equation. Thus, eleco
tric ﬁeld, physical conﬁnement of electrons, and in further macroscopic performances of
nano-MOSFETs will depend on these microscopic details. Phonon-electron interactions are
modeled as part of transport energy in a fashion of density functional theory, in which density
of phonons are assumed to be homogeneous in the channel and follows the Bose-Einstein statistics. Finally, the tunneling eﬀects of devices are measured as tunneling ratios of electron
currents.
Much of this dissertation has been written based on our papers:
• Langhua Hu and Guo-Wei Wei, Nonlinear Poisson equation for inhomogeneous media.
Biophysical Journal, 103, 758-766 (2012).
• Langhua Hu, Duan Chen and Guo-Wei Wei, High-order fractional partial diﬀerential
equations for molecular surface construction,Molecular Based Mathematical Biology,
1-25 (2012).
• Langhua Hu, Duan Chen and Guo-Wei Wei, Geometric and phonon eﬀects on nano
transistors, submitted (2013).
• Langhua Hu, Siyang Yang and Guo-Wei Wei, Simulation of inviscid compressible ﬂows
using PDE transform, submitted (2013).

166

5.2

Future work

In our proposed nonlinear Poisson equation (NLPE), we only consider the molecular charge
density, which is treated as the point partial charge source as showed in Eq. (1.9). It would
be an interesting to consider the salt eﬀect and incorporate the solvent charge distribution
in the nonlinear Poisson equation. When the solvent charge in the equilibrium, we can use
the Boltzmann distribution to approximate the solvent charge distribution, which is
Nc

ρs (r) =

ci qi e−φ(r)qi /kB T

(5.1)

i=1

where qi is the charge of ion species i, Nc is the number of ion species, kB is the Boltzmann
constant, T is the temperature, and ci is the bulk concentration of the ith ionic species.
When the solvent charge distribution given by Eq. (5.1) appears in the right hand side of
the nonlinear Poisson equation in Eq. (1.14), the NLPE will become much more nonlinear.
It remain an open question how the proposed nonlinear Poisson equation with salt eﬀect
will impact the electrostatic analysis. Furthermore, when the underlying chemical and biological system is far from the equilibrium, the Boltzmann distribution may no longer be
valid. Usually, these systems are described by coupled Nernst-Planck equation and Poisson equation, which is replaced by coupled Nernst-Planck equation and nonlinear Poisson
equation in the new framework. It’s a topic of interest to study how well the conjuncted
Nernst-Planck equation and nonlinear Poisson equation will work. Moreover, there are still
a huge room to explore in the ﬁeld of high order PDE transform and fractional PDE transform. For example, the molecular surface generation with the proposed high order fractional
PDE transform is much more time eﬃcient than MSMS approach, which makes it possible
to greatly facilitate the molecular dynamics. At last, it would be interesting to see how
167

the eﬀect of geometry,thermal and quantum will impact more complicated nano-transistors,
such as the nano-scale bio-sensor.

168

BIBLIOGRAPHY

169

BIBLIOGRAPHY
[1] A. Abrashkin, D. Andelman, and H. Orland. Dipolar poisson-boltzmann equation:
Ions and dipoles close to charge interfaces. Phys. Rev. Lett., 99:077801, 2007.
[2] O. P. Agrawal. Formulation of euler-lagrange equations for fractional variational problems. J. Math. Anal. Appl., 272:368 379, 2002.
[3] S. S. Ahmed, C. Ringhofer, and D. Vasileska. Parameter-free eﬀective potential method
for use in particle-based device simulations. IEEE Tran. Nanotech., 4:456 471, 2008.
[4] S. S. Ahmend and D. Vasileska. Modelling of narrow-width SOI devices. Semicond.
Sci. Technol, 19:S131–133, 2004.
[5] D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. Computational Fluid Mechanics
and Heat Transfer. McGraw-Hill, New York, 1984.
[6] P. Andrei and I. Mayergoyz. Analysis of ﬂuctuations in semiconductor devices through
self-consistent Poisson-Schr¨dinger computations. J. Applied Physics, 96:2071–2079,
o
2004.
[7] C. Azuara, H. Orland, M. Bon, P. Koehl, and M. Delarue. Incorporating dipolar
solvents with variable density in poisson-boltzmann electrostatics. Biophysical Journal,
95:5587C5605, 2008.
[8] B. Baeumer, M. Meerschaert, D. Benson, and S. Wheatcraft. Subordinated advectiondispersion equation for contaminant transport. Water Resour.Res., 37:1543–1550,
2001.
[9] J. Bai and X. C. Feng. Fractional-order anisotropic diﬀusion for image denoising. IEEE
Trans. Image Proc., 16:2492 2502, 2007.
[10] N. A. Baker. Poisson-Boltzmann methods for biomolecular electrostatics. Methods in
Enzymology, 383:94–118, 2004.
[11] N. A. Baker. Improving implicit solvent simulations: a Poisson-centric view. Current
Opinion in Structural Biology, 15(2):137–43, 2005.

170

[12] S. Barraud. Quantization eﬀects on the phonon-limited electron mobility in ultrathin
soi, ssoi and geoi devices. Semiconductor Scienec and Technology, 22:413–417, 2007.
[13] P. W. Bates, Z. Chen, Y. H. Sun, G. W. Wei, and S. Zhao. Geometric and potential
driving formation and evolution of biomolecular surfaces. Journal of Mathematical
Biology, 59(2):193–231, 2009.
[14] P. W. Bates, Z. Chen, Y. H. Sun, G. W. Wei, and S. Zhao. Geometric and potential
driving formation and evolution of biomolecular surfaces. J. Math. Biol., 59:193–231,
2009.
[15] P. W. Bates, G. W. Wei, and S. Zhao. The minimal molecular surface. arXiv:qbio/0610038v1, [q-bio.BM], 2006.
[16] P. W. Bates, G. W. Wei, and S. Zhao. Minimal molecular surfaces and their applications. Journal of Computational Chemistry, 29(3):380–91, 2008.
[17] S. Bednarek, B. Szafran, and J. Adamowski. Solution of the Poisson-Schr¨dinger
o
problem for a single-electron transistor. Zeitschrift fur Angewandte Mathematik und
Physik, 61:4461–4464, 2000.
[18] D. Beglov and B. Roux. Solvation of complex molecules in a polar liquid: an integral
equation theory. Journal of Chemical Physics, 104(21):8678–8689, 1996.
[19] A. L. Bertozzi and J. B. Greer. Low-curvature image simpliﬁers: Global regularity
of smooth solutions and laplacian limiting schemes. Communications on Pure and
Applied Mathematics, 57(6):764–790, 2004.
[20] F. Bianco, G. Puppo, and G. Russo. High-order central diﬀerencing for hyperbolic
conservation laws. SIAM J. Sci. Comput., pages 294–322, 1999.
[21] J. Blinn. A generalization of algebraic surface drawing. ACM Transactions on Graphics, 1(3):235–256, 1982.
[22] A. Blumen, G. Zumofen, and J. Klafter. Transport aspects in anomalous diﬀusion:
L’evy walks. Phys.Rev.A, 40:3964C3973, 1989.
[23] H. Borli, S. Kolberg, T. A. Fjeldly, and B. Iniguez. Precise modeling framework
for short-channel double-gate and gate-all-around mosfets. IEEE T Electron Devices,
55:2678–2686, 2008.

171

[24] I. Borukhov and D. Andelman. Steric eﬀects in electrolytes: A modiﬁed poissonboltzmann equation. Phys. Rev. Lett., 79(3):435–438, 1997.
[25] A. D. Buckingham and J. A. Pople. Theoretical studies of the kerr eﬀect i: Deviations
from a linear polarization law. Proc. Phys. Soc. A, 68:905– 909, 1955.
[26] F. M. Buﬂer, A. Schenk, and W. Fichtner. Eﬃcient monte carlo device modeling.
IEEE T Electron Devices, 47:1891–1897, 2000.
[27] W. Cai and C. W. Shu. Uniform high-order spectral methods for one- and twodimensional euler equations. J. Comput. Phys., 104:427–443, 1993.
[28] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods in ﬂuid
dynamics. Springer-Verlag, Berlin, 1988.
[29] M. Caputo. Linear model of dissipation whose w is almost frequency independent.
Geophys. J. R. Astr. Soc., 13:529 539, 1997.
[30] A. Chamberlin, C. Cramer, and D. Truhlar. Predicting aqueous free energies of solvation as functions of temperature. J. Phys. Chem. B, 110:5665–5675, 2006.
[31] A. Chambolle and P. L. Lions. Image recovery via total variation minimization and
related problems. Numerische Mathematik, 76(2):167–188, 1997.
[32] T. Chan, A. Marquina, and P. Mulet. High-order total variation-based image restoration. SIAM Journal on Scientiﬁc Computing, 22(2):503–516, 2000.
[33] D. Chen, Z. Chen, C. Chen, W. H. Geng, and G. W. Wei. MIBPB: A software package
for electrostatic analysis. J. Comput. Chem., 32:657 – 670, 2011.
[34] D. Chen, Z. Chen, and G. W. Wei. Quantum dynamics in continuum for proton transport II: Variational solvent-solute intersurface. International Journal for Numerical
Methods in Biomedical Engineering, DOI: 10.1002/cnm.1458, 2011.
[35] D. Chen and G. W. Wei. Modeling and simulation of electronic structure, material
interface and random doping in nano-electronic devices. J. Comput. Phys., 229:4431–
4460, 2010.
[36] F. Chen, C. M.and Liu, I. Turner, and V. Anh. A fourier method for the fractional
diﬀusion equation describing sub-diﬀusion. Journal of Computational Physics, 227:886
897, 2007.
172

[37] Z. Chen, N. A. Baker, and G. W. Wei. Diﬀerential geometry based solvation models
I: Eulerian formulation. J. Comput. Phys., 229:8231–8258, 2010.
[38] Z. Chen, N. A. Baker, and G. W. Wei. Diﬀerential geometry based solvation models
II: Lagrangian formulation. J. Math. Biol., 63:1139 1200, 2011.
[39] Z. Chen and G. W. Wei. Diﬀerential geometry based solvation models III: Quantum
formulation. J. Chem. Phys., 135(194108), 2011.
[40] Y. Cheng, I. M. Gamba, A. Majorana, and C.-W. Shu. A discontinuous galerkin
solver for boltzmannpoisson systems in nano devices. Computer Methods in Applied
Mechanics and Engineering, 198:3130–3150, 2009.
[41] Y. D. Cheng, I. M. Gamba, A. Majorana, and C. W. Shu. A discontinuous galerkin
solver for Boltzmann-Poisson systems in nano devices. Comput. Methods Appl. Mech.
Engrg., 198:3130–3150, 2009.
[42] S. Costolanski, Anne and C. Kelley. Eﬃcient Solution of the Wigner-Poisson Equation
for Modeling Resonant Tunneling Diodes. IEEE Transactions on Nanotechnology,
9:708–715, 2010.
[43] Y. Cui, Z. Zhong, D. Wang, W. U. Wang, and C. M. Lieber. High performance silicon
nanowire ﬁeld eﬀect transistors. NANO Letters, 3:149–152, 2003.
[44] S. Datta. Nanoscale device modeling: the Green’s function method. Superlattices and
Microstructures, 28:253–278, 2000.
[45] S. Datta. Electronic Transport in Mesoscopic Systems. Cambridge University Press,
1995.
[46] M. E. Davis and J. A. McCammon. Electrostatics in biomolecular structure and dynamics. Chemical Reviews, 94:509–21, 1990.
[47] C. de Falco, J. W. Jerome, and R. Sacco. A self-consitent iterative scheme for the onedimensional steady-state transistor calculations. IEEE Trans. Ele. Dev., 11:455–465,
1964.
[48] C. de Falco, J. W. Jerome, and R. Sacco. Quantum-corrected drift-diﬀusion models:
Solution ﬁxed point map and ﬁnite element approximation. J. Comput. Phys., 204:533–
561, 2009.

173

[49] S. Didas, J. Weickert, and B. Burgeth. Properties of higher order nonlinear diﬀusion
ﬁltering. Journal of mathematical imaging and vision, 35(3):208–226, 2009.
[50] T. J. Dolinsky, J. E. Nielsen, J. A. McCammon, and N. A. Baker. PDB2PQR: An
automated pipeline for the setup, execution, and analysis of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Research, 32:W665–W667, 2004.
[51] W. S. Don. Numerical study of pseudospectral methods in shock wave applications.
J. Comput. Phys, 110:103–111, 1994.
[52] W. S. Don and D. Gottlieb. Spectral simulation of supersonic reactive ﬂows. SIAM J.
Numer Anal., 35:2370–2384, 1998.
[53] F. Dong, B. Olsen, and N. A. Baker. Computational methods for biomolecular electrostatics. Methods in Cell Biology, 84:843–70, 2008.
[54] X. F. Duan, Y. Huang, Y. Ciu, J. F. Wang, and C. M. Lieber. Indium phosphide
nanowires as building blocks for nanoscale electronic and optoelectronics devices. Nature, 409:66–69, 2001.
[55] J. Dzubiella, J. M. J. Swanson, and J. A. McCammon. Coupling hydrophobicity,
dispersion, and electrostatics in continuum solvent models. Physical Review Letters,
96:087802, 2006.
[56] A. H. Elcock and J. A. McCammon. A continuum solvation model for studying protein
hydration thermodynamics at high temperatures. J. Phys. Chem. B, 101:9624–9634,
1997.
[57] H. Federer. Curvature Measures. Trans. Amer. Math. Soc., 93:418–491, 1959.
[58] M. Feig and C. L. Brooks III. Recent advances in the development and application of
implicit solvent models in biomolecule simulations. Curr Opin Struct Biol., 14:217 –
224, 2004.
[59] M. Ferrier, R. Clerc, G. Ghibaudo, F. Boeuf, and T. Skotnicki. Analytical model for
quantization on strained and unstrained bulk nMOSFET and its impact on quasiballistic current. Solid- State Electrons, 50:69–77, 2006.
[60] M. V. Fischetti. Master-equation approach to the study of electronic transport in small
semiconductor devices. Phys. Rev. B, 59:4901–4917, 1999.

174

[61] F. Fogolari, A. Brigo, and H. Molinari. The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology. Journal of Molecular Recognition,
15(6):377–92, 2002.
[62] B. Fornberg. A Practical Guide to Pseudospectral Methods. Cambridge University
Press, Cambridge, 1996.
[63] D. Frydel. Polarizable poisson-boltzmann equation: The study of polarizability eﬀects
on the structure of a double layer. J. Chem. Phys., 134:234704, 2011.
[64] R. Gabdoulline and R. Wade. Analytically deﬁned surfaces to analyze molecular interaction properties. Journal of Molecular Graphics, 14(6):341–353., 1996.
[65] P. Gaillardon, F. Clermidy, I. O’connor, J. Liu, M. Amadou, and G. Nicolescu. Matrix
nanodevice-based logic architectures and associated functional mapping method. ACM
Journal on Emerging Technologies in Computing Systems, 7:1–23, 2011.
[66] W. Geng and G. W. Wei. Multiscale molecular dynamics using the matched interface
and boundary method. J Comput. Phys., 230(2):435–457, 2011.
[67] W. H. Geng, S. N. Yu, and G. W. Wei. Treatment of charge singularities in implicit
solvent models. Journal of Chemical Physics, 127:114106, 2007.
[68] J. Giard and B. Macq. Molecular surface mesh generation by ﬁltering electron density
map. International Journal of Biomedical Imaging, 2010(923780):9 pages, 2010.
[69] G. Gilboa, N. Sochen, and Y. Y. Zeevi. Forward-and-backward diﬀusion processes for
adaptive image enhancement and denoising. IEEE Transactions on Image Processing,
11(7):689–703, 2002.
[70] G. Gilboa, N. Sochen, and Y. Y. Zeevi. Image sharpening by ﬂows based on triple well
potentials. Journal of Mathematical Imaging and Vision, 20(1-2):121–131, 2004.
[71] M. K. Gilson, M. E. Davis, B. A. Luty, and J. A. McCammon. Computation of electrostatic forces on solvated molecules using the Poisson-Boltzmann equation. Journal
of Physical Chemistry, 97(14):3591–3600, 1993.
[72] G. Giraud, C. M. Gordon, I. R. Dunkin, , and K. Wynne. The eﬀects of anion and
cation substitution on the ultrafast solvent dynamics of ionic liquids: A time-resolved
optical kerr-eﬀect spectroscopic study. J. Chem. Phys., 119:464–477, 2003.

175

[73] H. P. Gong and K. F. Freed. Langevin-debye model for nonlinear electrostatic screening
of solvated ions. Physical Review Letters, 102:057603, 2009.
[74] R. Gorenﬂo, F. Mainardi, E. Scalas, and M. Raberto. Fractional calculus and
continuous-time ﬁnance.iii,the diﬀusion limit.mathematical ﬁnance(konstanz, 2000).
Trends in Math., Birkhuser, Basel, page 171 C18, 2001.
[75] D. Gottlieb and J. S. Hesthaven. Spectral methods for hyperbolic problems. J. Comput.
Appl. Math., 128:83–131, 2001.
[76] D. Gottlieb, L. Lustman, and S. A. Orszag. Spectral calculations of one-dimensional
inviscid compressible ﬂow. SIAM J. Sci. Stat. Comput., 2:296–310, 1981.
[77] D. Gottlieb, C. W. Shu, A. Solomonoﬀ, and H. Vandeven. On the Gibbs phenomenon
I. J. Comput. Appl. Math., 43:81–98, 1992.
[78] D. Gottlieb and E. Tadmor. Recovering pointwise values of discontinuous data within
spectral accuracy. Progress and Supercomputing in Computational Fluid Dynamics,
Earll M. Murman and Saul S. Abarbanel eds.,Birkh¨user, pages 357–375, 1985.
a
[79] T. Grahs and T. Sonar. Entropy-controlled artiﬁcial anisotropic diﬀusion for the numerical solution of conservation laws based on algorithms from image processing. J.
Vis. Commun. Image R., 13:176–194, 2002.
[80] J. Grant and B. Pickup. A gaussian description of molecular shape. Journal of Physical
Chemistry, 99:3503–3510, 1995.
[81] J. B. Greer and A. L. Bertozzi. H-1 solutions of a class of fourth order nonlinear
equations for image processing. Discrete and Continuous Dynamical Systems, 10(12):349–366, 2004.
[82] J. B. Greer and A. L. Bertozzi. Traveling wave solutions of fourth order PDEs for
image processing. SIAM Journal on Mathematical Analysis, 36(1):38–68, 2004.
[83] P. Grochowski and J. Trylska. Continuum molecular electrostatics, salt eﬀects and
counterion binding. a review of the Poisson-Boltzmann theory and its modiﬁcations.
Biopolymers, 89(2):93–113, 2007.
[84] W. J. Gross, D. Vasileska, and D. K. Ferry. Three-dimensional simulations of ultrasmall
metal-oxide-semiconductor ﬁeld-eﬀect transistors: The role of the discrete impurities

176

on the device terminal characteristics. Journal of Applied Physics, 91(6):3737–3740,
2002.
[85] Y. Gu and G. W. Wei. Conjugated ﬁlter approach for shock capturing. Commun.
Numer. Meth. Engng., 19:99–110, 2003.
[86] Y. Gu, Y. C. Zhou, and G. W. Wei. Conjugate ﬁlters with spectral-like resolution for
2D incompressible ﬂows. Computers and Fluids, 33:777– 794, 2003.
[87] F. Hache, D. Ricard, C. Flytzanis, and U. Kreibig. The optical kerr eﬀect in small
metal particles and metal colloids: The case of gold. Appl. Phys. A, 47:437– 457, 1988.
[88] Z. Y. Han, N. Goldsman, and C. K. Lin. Incorporation of quantum corrections to
semiclassical two-dimensional device modeling with the Wigner-Boltzmann equation.
Solid-State Electronics, 49:145–154, 2005.
[89] U. Harbola and S. Mukamel. Superoperator nonequilibrium Green’s function theory of
many-body systems; applications to charge transfer and transport in open junctions.
Physics Report - Review Section of Physics Letters, 465:191–222, 2008.
[90] A. Harten, B. Engquist, S. Osher, and S. Chakravarthy. Uniform high-order accurate
essentially non-oscillatory schemes, III. J. Comput. Phys., 131:3–47, 1997.
[91] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:B864–871,
1964.
[92] B. Honig and A. Nicholls. Classical electrostatics in biology and chemistry. Science,
268(5214):1144–9, 1995.
[93] M. Y. Hussaini, D. A. Kopriva, M. D. Salas, and T. A. Zang. Spectral methods for
the euler equation: Part I - Fourier methods and shock-capturing. AIAA J., 23(No.
1):64–70, 1985.
[94] A. Jakalian, B. L. Bush, D. B. Jack, and C. I. Bayly. Fast, eﬃcient generation of
high-quality atomic charges. am1-bcc model: I. method. Journal of Computational
Chemistry, 21(2):132–146, 2000.
[95] G. S. Jiang and C. W. Shu. Eﬃcient implementation of weighted ENO schemes. J.
Comp. Phys., 126(1):202–228, 1996.

177

[96] H. Y. Jiang, S. H. Shao, W. Cai, and P. W. Zhang. Boundary treatments in
non-equilibrium Green’s function (NEGF) methods for quantum transport in nanoMOSFETs. J Comput. Physics, 227:6553–6573, 2008.
[97] X. Jiang, H. X. Deng, J. W. Luo, S. S. Li, and L. W. Wang. A fully three-dimensional
atomistic quantum mechanical study on random dopant-induced eﬀects in 25-nm MOSFETs. IEEE T Electron Devices, 55:1720–1726, 2008.
[98] S. Jin and C. D. Levermore. Numerical schemes for hyperbolic conservation laws with
stiﬀ relaxation terms. J. Comput. Phys., 126:449–467, 1996.
[99] Z. M. Jin and X. P. Yang. Strong solutions for the generalized Perona-Malik equation for image restoration. Nonlinear Analysis-Theory Methods and Applications,
73(4):1077–1084, 2010.
[100] B. Jogai. Inﬂuence of surface states on the two-dimensional electron gas in algan/gan
heterojunction ﬁeld-eﬀect transistors. J. Applied Physics, 93:1631–1635, 2003.
[101] S. T. Johansen, J. Y. Wu, and W. Shyy. Filter-based unsteady rans computations.
Int. J. Heat Fluid Flow, 25:10–21, 2004.
[102] M. Kanduc, M. A. Naji, Y. S. Jho, P. A. Pincus, and R. Podgornik. The role of
multipoles in counterion-mediated interactions between charged surfaces: strong and
weak coupling. Journal of Physics: Condensed Matter, 21:424103, 2009.
[103] H. R. Khan, D. Vasileska, and S. S. Ahmed. Modeling of FinFet: 3D MC Simulation Using FMM and Unintentional Doping Eﬀects on Device Operation. Journal of
Computational Electronics, 3:337–340, 2004.
[104] T. Khan, D. Vasileska, and T. J. Thornton. Eﬀect of interface roughness on silicon-oninsulator-metal-semiconductor ﬁeld-eﬀect transistor mobility and the device low-power
high-frequence operation. J. Vac. Sci. Technol. B, 23(4):1782–1784, 2005.
[105] J. G. Kirkwood. Theory of solution of molecules containing widely separated charges
with special application to zwitterions. J. Comput. Phys., 7:351 – 361, 1934.
[106] P. Koehl, H. Orland, and M. Delarue. Beyond the poisson-boltzmann model: Modeling
biomolecule-water and water-water interactions. Phys. Rev. Lett., 102:087801, 2009.
[107] W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation
eﬀects. Phys. Rev., 140:A1133– A1138, 1965.
178

[108] A. Kolpakov, A. K. Tagantsev, L. Berlyand, and A. Kanareykin. Nonlinear dielectric
response of periodic composite materials. J Electroceram, 18:129–137, 2007.
[109] A. A. Kornyshev and G. Sutmann. Nonlocal dielectric saturation in liquid water.
Physical Review Letters, 79:3435–3438, 1997.
[110] R. Krasny. A study of singularity formation in a vortex sheet by the point vortex
approximation. J. Fluid Mech., 167:65–93, 1986.
[111] A. Kurganov and D. Levy. A third-order semidiscrete central scheme for conservation
laws and convention diﬀusion equations. SIAM J. Sci. Comput., 22:1461–1488, 2000.
[112] S. Lee, S. K. Lele, and P. Moin. Interaction of isotropic turbulence with shock waves:
Eﬀect of shock strength. J. Fluid Mech., 340:225–247, 1997.
[113] R. J. LeVeque and M. Pelanti. A class of approximate riemann solvers and their
relation to relaxation schemes. J. Comput. Phys., 172:574–591, 2001.
[114] R. M. Levy, L. Y. Zhang, E. Gallicchio, and A. K. Felts. On the nonpolar hydration free
energy of proteins: surface area and continuum solvent models for the solute-solvent
interaction energy. Journal of the American Chemical Society, 125(31):9523–9530,
2003.
[115] D. R. Lide. CRC Handbook of Chemistry and Physics, volume 75. CRC:Boca Raton,
FL, 1995.
[116] X. D. Liu and E. Tadmor. Third order non-oscillatory central scheme for hyperbolic
conservation laws. Numer. Math., 79:397–425, 1998.
[117] B. J. Loughnane, A. Scodinu, R. A. Farrer, J. T. Fourkas, and U. Mohanty. Exponential
intermolecular dynamics in optical kerr eﬀect spectroscopy of small-molecule liquids.
J. Chem. Phys., 111:2686–2694, 1999.
[118] J. H. Luscombe, A. M. Bouchard, and M. Luban. Electron conﬁnement in quantum nanostructure - Self-consistent Poisson-Schr¨dinger theory. Physical Review B,
o
46:10262–10268, 1992.
[119] M. Lysaker, A. Lundervold, and X. C. Tai. Noise removal using fourth-order partial
diﬀerential equation with application to medical magnetic resonance images in space
and time. IEEE Transactions on Image Processing, 12(12):1579–1590, 2003.

179

[120] J. MacKerell, A. D., D. Bashford, M. Bellot, J. Dunbrack, R. L., J. D. Evanseck, M. J.
Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera,
F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, I. Reiher,
W. E., B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe,
J. Wiorkiewicz-Kuczera, D. Yin, and M. Karplus. All-atom empirical potential for
molecular modeling and dynamics studies of proteins. Journal of Physical Chemistry
B, 102(18):3586–3616, 1998.
[121] F. Mainardi and R. Gorenﬂo. On Mittag-Leﬄer-type functions in fractional evolution
processes. Journal of Computational and Applied Mathematics, 118:283 – 299, 2000.
[122] M. H. Majles Ara, S. H. Mousavi, Z. Mousavi, and M. S. Zakerhamidi. Investigation
of the kerr eﬀect and third-order susceptibility constants in a nematic liquid crystal.
Journal of Molecular Liquids, 161:41 – 43, 2011.
[123] J. F. McKenzie and K. O. Westphal. Interaction of linear waves with oblique shock
waves. Physics of Fluids, 11:2350–2362, 1968.
[124] W. R. McKinnon. Fokker-Planck approach to extending the one-ﬂux method of carrier
transport in semiconductors to variable energies. J. Applied Phys., 94:4986–4994, 2003.
[125] M. Meerschaert. Fractional calculus, anomalous diﬀusion, and probability. Fractional
Dynamics, R. Metzler and J. Klafter, Eds., World Scientiﬁc, Singapore, pages 265–
284, 2012.
[126] M. Meerschaert and C. Tadjeran. Finite diﬀerence approximations for fractional
advection-dispersion ﬂow equations. Journal of Computational and Applied Mathematics, 172(1):65–77, 2004.
[127] B. Mennucci, E. Cances, and J. Tomasi. Evaluation of solvent eﬀects in isotropic and
anisotropic dielectrics and in ionic solutions with a uniﬁed integral equation method:
theoretical bases, computational implementation, and numerical applications. J. Phys.
Chem. B, 101:10506 – 10517, 1997.
[128] H. Nessyahu and E. Tadmor. Non-oscillation central diﬀerencing for hyperbolic conservation laws. J. Comput. Phys., 87:408–463, 1990.
[129] A. Nicholls, D. L. Mobley, P. J. Guthrie, J. D. Chodera, and V. S. Pande. Predicting small-molecule solvation free energies: An informal blind test for computational
chemistry. Journal of Medicinal Chemistry, 51(4):769–79, 2008.

180

[130] M. Nina, D. Beglov, and B. Roux. Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulations. J. Phys. Chem. B,
101:5239–5248, 1997.
[131] P. Nithiarasu, O. C. Zienkiewicz, B. V. K. S. Sai, K. Morgan, R. Codina, and
M. Vazquez. Shock capturing viscosities for the general ﬂuid mechanics algorithm.
Int. J. Numer. Meth. Fluids, 28:1325–1353, 1998.
[132] A. Papazyan and A. Warshel. A stringent test of the cavity concept in continuum
dielectrics. J. Chem. Phys., 107:7975– 7978, 1997.
[133] R. Parr and W. Yang. Density Functional Theory of Atoms and Molecules. 1989.
[134] P. Perona and J. Malik. Scale-space and edge-detection using anisotropic diﬀusion.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629–639, 1990.
[135] J. X. Qiu and C. W. Shu. On the construction, comparison, and local characteristic
decomposition for high order central WENO schemes. J. Comput. Phys., 183:187–209,
2002.
[136] M. Raberto, E. Scalas, and F. Mainardi. Waiting-times and returns in high-frequency
ﬁnancial data: an empirical study. Physica A, 314:749C755, 2002.
[137] J. Robertson and D. Hall. Nonlinear dielectric properties of particulate barium
titanate-polymer composites. Journal of Physics D: Applied Physices, 41:115407–
115414, 2008.
[138] W. Rocchia, E. Alexov, and B. Honig. Extending the applicability of the nonlinear
poisson-boltzmann equation: Multiple dielectric constants and multivalent ions. J.
Phys. Chem., 105:6507–6514, 2001.
[139] B. Roux and T. Simonson. Implicit solvent models. Biophysical Chemistry, 78(1-2):1–
20, 1999.
[140] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal
algorithms.
[141] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. In Proceedings of the eleventh annual international conference of the Center
for Nonlinear Studies on Experimental mathematics : computational issues in nonlinear science, pages 259–268, Amsterdam, The Netherlands, The Netherlands, 1992.
Elsevier North-Holland, Inc.
181

[142] L. Sabatelli, S. Keating, J. Dudley, and P. Richmond. Waiting time distributions in
ﬁnancial markets. Eur.Phys.J.B, 27:273C275, 2002.
[143] M. F. Sanner, A. J. Olson, and J. C. Spehner. Reduced surface: An eﬃcient way to
compute molecular surfaces. Biopolymers, 38:305–320, 1996.
[144] S. Scaldaferri, G. Curatola, and G. Iannaccone. Direct solution of the Boltzmann
transport equation and Poisson-Schr¨dinger equation for nanoscale MOSFETs . IEEE
o
T. Electron Devices, 54:2901–2909, 2007.
[145] M. J. Schnieders, N. A. Baker, P. Ren, and J. W. Ponder. Polarizable atomic multipole
solutes in a Poisson-Boltzmann continuum. Journal of Chemical Physics, 126:124114,
2007.
[146] K. A. Sharp and B. Honig. Calculating total electrostatic energies with the nonlinear
Poisson-Boltzmann equatlon. Journal of Physical Chemistry, 94:7684–7692, 1990.
[147] T. Simonson. Macromolecular electrostatics: continuum models and their growing
pains. Current Opinion in Structural Biology, 11(2):243–252, 2001.
[148] T. Simonson. Electrostatics and dynamics of proteins. Reports on Progress in Physics,
66(5):737–87, 2003.
[149] J. C. Slater. A simpliﬁcation of the Hartree-Fock method. Phys. Rev., 81:385–390,
1951.
[150] R. F. Snider. Quantum-mechanical modiﬁed boltzmann equation for degenerate internal states. J. Chem. Phys., 32:1051–1060, 1960.
[151] Y. H. Sun, P. R. Wu, G. W. Wei, and G. Wang. Evolution-operator-based single-step
method for image processing. Int. J. Biomed. Imaging, 83847:1–27, 2006.
[152] Y. H. Sun, Y. C. Zhou, S. G. Li, and G. W. Wei. A windowed fourier pseudospectral
method for hyperbolic conservation laws. J. Comput. Phys., 214:466–490, 2006.
[153] J. M. J. Swanson, J. Mongan, and J. A. McCammon. Limitations of atom-centered
dielectric functions in implicit solvent models. Journal of Physical Chemistry B,
109(31):14769–72, 2005.
[154] J. M. J. Swanson, J. A. Wagoner, N. A. Baker, and J. A. McCammon. Optimizing the
poisson dielectric bounday with explicit solvent forces and energies: Lessons learned
182

with atom-centered dielectric functions. Journal of Chemical Theory and Computation,
3(1):170–83, 2007.
[155] E. Tadmor. Convergence of spectral methods for nonlinear conservation laws. SIAM
J. Numer. Anal., 26:30–44, 1989.
[156] A. K. Tagantsev, V. O. Sherman, K. F. Astaﬁev, J. Venkatesh, and N. Setter. Ferroelectric materials for microwave tunable applications. Journal of Electroceramics,
11:5–66, 2003.
[157] T. Tasdizen, R. Whitaker, P. Burchard, and S. Osher. Geometric surface processing
via normal maps. Acm Transactions on Graphics, 22(4):1012–1033, 2003.
[158] E. F. Toro, R. C. Millington, and V. A. Titarev. Der: arbitrary-order non-oscillatory
advection scheme. Proc. 8th International Conference on Non-linear Hyperbolic Problems, Magdeburg, Germany, March, 2000.
[159] L. N. Trefethen. Spectral Methods in Matlab. Oxford University, Oxford, England,
2000.
[160] A. Trellakis, A. T. Galick, A. Pacelli, and U. Ravaioli. Iteration scheme for the solution
of the two-dimensional Schr¨dinger-Poisson equations in quantum structures. J. Appl.
o
Phys., 81:7880–7884, 1997.
[161] V. C. Tung, M. J. Allen, Y. Yang, and R. B. Kaner. High-throughput solution processing of large-scale graphene. Nature Nanotechnology, 4:25–29, 2009.
[162] H. Vandeven. Family of spectral ﬁlters for discontinuous problems. J. Sci. Comput.,
6:159–192, 1991.
[163] P. Vath, M. B. Zimmt, D. V. Matyushov, and G. A. Voth. A failure of continuum theory: Temperature dependence of the solvent reorganization energy of electron transfer
in highly polar solvents. J. Phys. Chem. B, 103:9130– 9140, 1999.
[164] U. V.N. and P. E.A., editors. Methods in Protein Structure and Stability Analysis:
Conformational Stability, Size, Shape and Surface of Protein Molecules. Nova sicence
publishers inc., 2007.
[165] J. A. Wagoner and N. A. Baker. Assessing implicit models for nonpolar mean solvation
forces: the importance of dispersion and volume terms. Proceedings of the National
Academy of Sciences of the United States of America, 103(22):8331–6, 2006.
183

[166] L. Waldmann. Die Boltzmann-Gleichung fur Gase mit rotierenden Molekulen. Z.
Naturforsch. Teil A, 12:660–662, 1957.
[167] D. C. Wan, B. S. V. Patnaik, and G. W. Wei. Discrete singular convolution-ﬁnite
subdomain method for the solution of incompressible viscous ﬂows. J. Comput. Phys.,
180:229–255, 2002.
[168] J. Wang, E. Polizzi, and M. Lundstrom. A three-dimensional quantum simulation of
silicon nanowire transistors with the eﬀective-mass approximation. J. Appl. Phys.,
96:2192–203, 2004.
[169] Y. Wang, G. W. Wei, and S.-Y. Yang. Partial diﬀerential equation transform – Variational formulation and Fourier analysis. International Journal for Numerical Methods
in Biomedical Engineering, 27:1996–2020, 2011.
[170] Y. Wang, G. W. Wei, and S.-Y. Yang. Iterative ﬁltering decomposition based on local
spectral evolution kernel. Journal of Scientiﬁc Computing, pages DOI: 10.1007/s10915–
011–9496–0, accepted, 2011.
[171] Y. Wang, G. W. Wei, and S.-Y. Yang. Mode decomposition evolution equations.
Journal of Scientiﬁc Computing, accepted,2011.
[172] A. Warshel and L. ˚qvist. Electrostatic energy and macromolecular function. Annu.
A
Rev. Biophys. Biophys. Chem., 20:267–98, 1991.
[173] A. Warshel and A. Papazyan. Electrostatic eﬀects in macromolecules: fundamental
concepts and practical modeling. Current Opinion in Structural Biology, 8(2):211–7,
1998.
[174] A. Warshel, P. K. Sharma, M. Kato, and W. W. Parson. Modeling electrostatic effects in proteins. Biochimica et Biophysica Acta (BBA) - Proteins & Proteomics,
1764(11):1647–76, 2006.
[175] J. Warwicker and H. C. Watson. Calculation of the electric potential in the active site
cleft due to alpha-helix dipoles. Journal of Molecular Biology, 157(4):671–9, 1982.
[176] J. D. Weeks, D. Chandler, and H. C. Andersen. Role of repulsive forces in determining
the equilibrium structure of simple liquids. Journal of Chemical Physics, 54(12):5237–
47, 1971.
[177] G. W. Wei. Discrete singular convolution for the solution of the Fokker-Planck equations. J. Chem. Phys., 110:8930–8942, 1999.
184

[178] G. W. Wei. Generalized Perona-Malik equation for image restoration. IEEE Signal
Processing Letters, 6(7):165–167, 1999.
[179] G. W. Wei. Oscillation reduction by anisotropic diﬀusions. Comput. Phys. Commun.,
144:317–342, 2002.
[180] G. W. Wei. Diﬀerential geometry based multiscale models. Bulletin of Mathematical
Biology, 72:1562 – 1622, 2010.
[181] G. W. Wei and Y. Gu. Conjugated ﬁlter approach for solving Burgers’ equation. J.
Comput. Appl. Math., 149:439–456, 2002.
[182] G. W. Wei and Y. Q. Jia. Synchronization-based image edge detection. Europhysics
Letters, 59(6):814–819, 2002.
[183] G. W. Wei, Y. B. Zhao, and Y. Xiang. Discrete singular convolution and its application
to the analysis of plates with internal supports, I theory and algorithm. Int. J. Numer.
Meth. Engng., 55:913–946, 2002.
[184] G. W. Wei, Q. Zheng, Z. Chen, and K. Xia. Diﬀerential geometry based ion transport
models. SIAM Review, accepted, 2012.
[185] T. P. Witelski and M. Bowen. ADI schemes for higher-order nonlinear diﬀusion equations. Applied Numerical Mathematics, 45(2-3):331–351, 2003.
[186] A. Witkin. Scale-space ﬁltering: A new approach to multi-scale description. Proceedings of IEEE International Conference on Acoustic Speech Signal Processing, 9:150–
153, 1984.
[187] K. Xu. A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection
with artiﬁcial dissipation and Godunov method. J. Comput. Phys., 171:289–335, 2001.
[188] M. Xu and S. L. Zhou. Existence and uniqueness of weak solutions for a fourth-order
nonlinear parabolic equation. Journal of Mathematical Analysis and Applications,
325(1):636–654, 2007.
[189] Y. Xu, K. Zhu, S. Yan, Z. Jin, Y. Wang, H. Chen, J. Luo, and B. Yu. Quantum and
thermo-mechanical noise squeezing in nanoresonators: A comparative study. Applied
Physics Letters, 100:023105, 2012.

185

[190] H. C. Yee, N. D. Sandham, and M. J. Djomehri. Low-dissipative high-order shockcapturing methods using characteristic-based ﬁlters. J. Comput. Phys., 150:199–238,
1999.
[191] Y. You and M. Kaveh. Fourth-order partial diﬀerential equations for noise removal.
IEEE Transactions on Image Processing, 9(10):1723–1730, 2002.
[192] B. Yu, L. Chang, S. Ahmed, H. H. Wang, S. Bell, C. Yang, C. Tabery, C. Ho, Q. Xiang,
T. King, J. Bokor, C. Hu, M. Lin, and D. Kyser. Finfet scaling to 10 nm gate length.
IEDM Tech. Dig., pages 251–254, 2002.
[193] B. Yu, C. H. J. Wann, E. D. Nowak, K. Noda, and C. M. Hu. Short-channel effect improved by lateral channel-engineering in deep-submicronmeter mosfet’s. IEEE
TRANSACTIONS ON ELECTRON DEVICES, 44:627–634, 1997.
[194] S. N. Yu, W. H. Geng, and G. W. Wei. Treatment of geometric singularities in implicit
solvent models. Journal of Chemical Physics, 126:244108, 2007.
[195] S. N. Yu and G. W. Wei. Three-dimensional matched interface and boundary (MIB)
method for treating geometric singularities. J. Comput. Phys., 227:602–632, 2007.
[196] S. N. Yu, Y. C. Zhou, and G. W. Wei. Matched interface and boundary (MIB) method
for elliptic problems with sharp-edged interfaces. J. Comput. Phys., 224(2):729–756,
2007.
[197] G. Zaslavsky. Fractional kinetic equation for hamiltonian chaos.chaotic advection,
tracer dynamics and turbulent dispersion. Phys.D, 76:110C122, 1994.
[198] Y. Zhang, G. Xu, and C. Bajaj. Quality meshing of implicit solvation models of
biomolecular structures. Computer Aided Geometric Design, 23(6):510–30, 2006.
[199] S. Zhao and G. W. Wei. High-order FDTD methods via derivative matching for
Maxwell’s equations with material interfaces. J. Comput. Phys., 200(1):60–103, 2004.
[200] S. L. Zhao, R. Ramirez, R. Vuilleumier, and D. Borgis. Molecular density functional theory of solvation: from polar solvents to water. Journal of chemical physics,
134(19):194102, 2011.
[201] Q. Zheng, D. Chen, and G. W. Wei. Second-order Poisson-Nernst-Planck solver for
ion transport. Journal of Comput. Phys., 230:5239 – 5262, 2011.

186

[202] Q. Zheng and G. W. Wei. Poisson-Boltzmann-Nernst-Planck model. Journal of Chemical Physics, 134:194101, 2011.
[203] Q. Zheng, S. Yang, and G.-W. Wei. Molecular surface generation using pde transform.international journal for numerical methods in biomedical engineering. in press,
2011.
[204] Q. Zheng, S. Y. Yang, and G. W. Wei. Molecular surface generation using pde transform. International Journal for Numerical Methods in Biomedical Engineering, 2012.
[205] Y. C. Zhou, Y. Gu, and G. W. Wei. Shock-capturing with natural high frequency
oscillations. Int. J. Numer. Methods Fluid, 41:1319–1338, 2003.
[206] Y. C. Zhou and G. W. Wei. High resolution conjugate ﬁlters for the simulation of
ﬂows. J. Comput. Phys., 189:159–179, 2003.
[207] Y. C. Zhou and G. W. Wei. On the ﬁctitious-domain and interpolation formulations
of the matched interface and boundary (MIB) method. J. Comput. Phys., 219(1):228–
246, 2006.
[208] Y. C. Zhou, S. Zhao, M. Feig, and G. W. Wei. High order matched interface and
boundary method for elliptic equations with discontinuous coeﬃcients and singular
sources. J. Comput. Phys., 213(1):1–30, 2006.

187