Transcript 並列コンピューティングを援用した強結合分子動力学シミュレーターの高速

Tight-binding molecular dynamics study
of mechanical and electronic properties
in twisted graphene nanoribbons
Satofumi Souma, Shozo Kaino, and Matsuto Ogawa
Department of Electrical and Electronic Engineering,
Kobe University, Japan
SISPAD 2012 Denber, Colorado, USA, Sep. 6, 2012
Outline
Introduction
Calculation method
Numerical results
Summary
1
Introduction
•Recent rapid progress in nanostructure electronics based on
carbon based nano-scale materials, such as CNT, Graphene.
--- electronic property is sensitive to the geometrical structure.
It is essentially important to know the relaxed geometrical
structure for reliable simulations.
#First principles calculation (accurate, but time-consuming),
#Tight-binding molecular dynamics method.
#Classical valence force field method (Tersoff, Keating, etc..)
(suitable for large-scale simulation)
Efficient algorithm is required to
perform the structural relaxation
calculation for large scale system.
2
Purpose of this study
Graphene:
High electron mobility and thermal conductivity
However– No bandgap at Fermi energy
Essentially required to open the
bandgap in graphene for switching
device application.
example・・・
# Graphene nanoribbon(GNR) (especially AGNR)
Our question:
#Applying electric field to bi-layer graphene
#In-plane/out-of-plane strain
Can we further control the electronic state of AGNR (especially bandgap) by
applying various types of strain? (Stretching/Twisting)
Electronic structure is sensitive to the atomic size geometry
Essentially important to have the relaxed geometry under the strain
We use Tight-Binding Molecular Dynamics (TBMD) ） C. H. Xu, et al, J. Phys.
Condens. Matter . 4, 28 (1992).
3
Bandgap in AGNR
x
y
z
In this study, we choose 7×22, 9×22, 11×22 AGNR to analyze the
effect of streatching/twisting
7x22
9x22
11x22
4
Flowchart of structural relaxation calc.
based on TBMD
Setup the TB Hamiltonian for initial atomic positions
Total energy (band structure energy + repulsion energy )
Etot(=Ebs+Erep)
Force Ftot(tn) acting on atoms
Atomic positions ri(tn+1) are updated by using velocity Verlet method.
Force Fi(tn+1) and velocity vi(tn+1) are updated
Velocity Verlet method
5
Conventional diagonalization method
In TBMD method, bandstructure energy Ebs is a
sum of occupied energy eigenvalue for H
In finite temperature
In conventional method, εi iis calculated by
diagonalizing the Hamiltonian matrix H
However….
Calclation time of εi is proportional to the cupic of N (number
of atoms)  time-consuming
Inappropriate for parallel calculation.
We employ the Fermi operator expansion
method for the calculation of Ebs
6
Fermi Operator Expansion method
In FOE method, the bandstructure energy Ebs is calculated as
requires the matrix multiplication of H (Hamiltonian) and Fermi operator
Fermi distribution func.
In FOE, the Fermi distribution function is expressed by Chebyshev polynomial
expansion p(H), which enables us to obtain the matrix form of Fermi func.
7
Chebyshev polynomial
 related to the computational accuracy
In Chebyshev expansion, the polynomial p(H) satisfies the
following equation
Matrix multiprication H×Tｊ(H)
n
 H11
H
 21
T j 1   H 31

 
 H n1
H12
H
H13  H1n 






n
Tj
T11j T12j T13j  T 1nj  
 j
 
T21
 
T31j
 n

 

 
T j
 
 n1

T j 1
T11j 1 T12j 1 T13j 1  T j1n1 
 j 1

T21


 T31j 1


 

T j 1

 n1




n



nn
For 1 column of Tｊ+1(H) n×n, which has to be repeated for n-rows = n3
By using the sparse nature of H + cut-off radius for matrix
element, n3 is reduced to n.
8
Comparison between diagonalization and FOE
Typical example:
AGNR with 14 atoms in 1 unit cell.
Number of unit cell is increased.
Estimation of the time required to obtain
the relaxed structure.
Significant reduction
of calculation time
9
Result 1: Stretched AGNR
Fracture limit
SR dependence of total energy
Stretching ratio
SR 
L' L
100[%]
L
SR dendence of the bandgap
Stretched AGNR
←SR dependence of the bandgap
in stretched AGNR
Energy dependence of the
transmission for various SR
Transport gap and electrode
bandgap can be different
Electrode band structure
Result 2: Twisted AGNR
L
Twisting angle θ is assumed
to be controlled while the
ribbon length L is fixed.
Twisting causes the stretching especially near
the edge of ribbon
It is interesting to compare twisting and
stretching to understand the effect of
twisting on the property of AGNR.
Twisting angle θ 
Stretching ratio: SR ?
L’
SR 
L' L
100[%]
L
12
Fracture limit of twisted AGNR
How does the twisted AGNR look like?
ParameterτW
τW： quantifies how much the AGNR
twisted along the width direction with the
curvature along the length direction.
・θ[rad] ⇔ τW ⇔ SR[%]
Pekka Koskinen, 2011, April. Phys. Lett. 99, 013105
This relationship is determined so that the strain energy is comparable in
stretched and twisted cases
a. 9×22 AGNR 900°(τW =1.62), Tube-like
b. 7×22 AGNR 900°(τW =1.21), Ribbon-like
c. 7×22 AGNR 1210°(τW =1.62) , Tube-like
Twisted ribbons with the same τW look similarly
15
Fracture limit of twisted AGNR
 max
Fracture limit is roughly
proportional to the aspect ratio
L

W
GNR can be twisted more when
the length L is much longer than the
width W (consistent with intuition).
 max W  

W W
L
 max  
L
W
 max
W
W L
 max W 
W   max    
L
L
L W
Fracture limit (τW) is constant
(~2.6 [rad]) for L/W>6 (L>6W)
16
Twisting and Stretching 1
At τW=1.55～1.70 twisted AGNR changes its
relaxed geometry from ribbon like to tube like,
by reducing the total energy.
17
Twisting and Stretching 2
Twisting and twisting agree well until τW＝1.55
, and after that they start to deviate
Geometrical structure changes from ribbon-like
to tube-like, and then the large curvature
influences the electronic band structure.
7×22 AGNR
9×22 AGNR
11×22 AGNR
リ
ボ
ン
・ stretching
・ twisting
τW=1.55
τW=1.55
τW=1.62
チ
ュ
ー
ブ
的
τW=1.62
τW=1.62
τW=1.70
18
Summary
Summary
We employed order-N FOE method for TBMD structural relaxation
calculation along with the parallel calculation, and succeeded to
reduce the computation time significantly.
In Twisted AGNR, the fracture limit is determined by the aspect ratio.
Bandgap of AGNR can be controlled by stretching and twisting.
Twisted AGNR changes its relaxed structure from ribbon-like to
tube-like abruptly at around a critical curvature, then the electronic
property also changes.
Future plan
More care study on the effect of ribbon-like-to-tube-like transition
to the electronic property.
Similar study for zigzag edged GNR.
19
damping factor
• 多項式近似を行うと振動がみられるが、この振動をどのようにして抑えているか？
⇒ Chebyshev多項式近似にJackson damping factorを導入。これによりTBMDでより少
ないtime stepで収束するようになった。
1,2
Chebyshev
1
Chebyshev-Jackson
f(E,T)
0,8
0,6
M=75
0,4
0,2
0
-0,2
-10
-5
0
Energy (eV)
5
10
20
対角化法-FOEの誤差
⇒mpl=200として計算
C2分子について、最近接原子間距離a0を変化させたときの
全エネルギーEtotを計算
a0=1.413Åのときの結合エネルギーEbsを比較
対角化法
FOE O(N)
-22.113 eV
-22.111 eV
誤差： 0.002 eV/atom
誤差率： 0.009%
21
Ebandの計算
温度による影響を
取り入れると
5
n5=0
4
n4=0
3
n3=1
EF
2
n2=1
i=1
n1=1
22
Etotの計算
本研究では原子の持つ全エネルギーEtotを以下の形で与える
Etot=Ebs+Erep
Ebs： 結合エネルギー, 電子によって占有されているエネルギー固有値を足し合わせたもの
Erep： 原子が近づくことで生じる反発力を生じさせるエネルギー
Ebs, Erepとも2原子間の距離に依存する
C2分子
各エネルギーの関係
23
TBMDシミュレーション
(4,0)CNT-2layer 32atom
外力を加えてない状態の(4,0)CNTに対して
構造緩和を適用
Total Energy [eV/atom]
-6,9
-6,94
-6,98
全エネルギーEtotが最も低くな
るような安定構造に収束する
-7,02
Etot before： -6.910828 [eV/atom]
Etot after ： -7.079077 [eV/atom]
-7,06
-7,1
0
100
200
300
400 500 600
Time Step[ps]
700
800
900 1000
24
構造緩和にかかる時間
(4,0)CNT TBMD
4000
CPU time (s)
3500
y = 0,0098x2,6099
3000
2500
構造緩和にかかる時間は原子数N
の3乗に比例して大きくなる
2000
1500
1000
500
0
0
20
40
60
80
number of atoms
100
120
140
Ebsの計算にかかる時間が原子数N
の3乗に比例しており、計算量が肥大
してしまう原因となっている
(4,0)CNT Ebs
1200
y = 0,0002x3,21
CPU time (s)
1000
Ebsの計算時間を減らす
800
600
400
200
0
0
20
40
60
80
number of atoms
100
120
140
25
オーダーN法
オーダーN法とは・・・
通常のバンド計算手法では、扱う系の原子の数Nの3乗に比例して計算量が増えるが
これをNの１乗のオーダー（オーダーN）にしようとする電子状態計算手法
本研究ではオーダーN法である
Fermi operator expansion(FOE)を取り入れた
FOEのメリット
計算量が原子数Nの1乗に比例するため、対角化法と比較して計
算量を大幅に減らすことができ、巨大な系を扱うことが出来る。
並列計算機で容易に実行できる
26
フェルミ行列の多項式近似
ε→Hに置き換えると
近似の精度を
決定する
27
Step1. ハミルトニアンのスパース性
TBハミルトニアン行列Hは各要素にScaling関数s(rij)を組み込ん
でいるため、行列要素の値が原子間距離に依存している
行列要素のほとんどが0であるスパース行列（疎行列）になっている
n off
T j 1 
H
1 2 3 ・・・

1
① ② 

2


3 

 

・
・

・
n
Tj
T11j T12j T13j  T 1nj  
 j
 
T21
 
T31j
 n

 

 
T j
 
 n1

行列要素が原子番号に対応
T j 1
T11j 1 T12j 1 T13j 1  T j1n1 
 j 1

T21

j 1


 T31


 

T j 1

 n1




n



①
1
原子軌道
3
②
1
5
n off  n
スパース性を考慮して非ゼロの成分だけ計算に用いることで
計算量を減らせる（noff× n2, noff は原子数によらない一定の数）
28
Step2. cutoff半径rlocの導入
n off
1 2 3
1

2
3 
T j 1 
n



n
Tj
H
1







T j 1
2 3
T T12j T13j  T 1nj 


T

T





T j

 n1








j
11
j
21
j
31







n loc
n off  n loc
行列要素間の物理的距離に対してcutoffを設ける
上の図の場合、a1-n＞ rlocならば（Hのn列）×(Tjの1列)の計算を丸々無視できる
⇒Tj+1の要素をnloc個だけ求める
rloc
a1-n
1
n
nlocは原子数によらない一定の数なので、
ハミルトニアンのスパースと組み合わせると全計算量はnoff× nloc×nでオーダーNとなる
29