計算数理工学論文集 Vol. 7, No.2 (2008 年 3 月), 論文 No. 18-080317

JASCOME

内部自由度を有する結晶のための第一原理密度汎関数法に基づく
高速局所擬似連続体法
DATABASE-ACCELERATED PARALLELIZED LOCAL QUASICONTINUUM METHODBASED ON FIRST-PRINCIPLES DENSITY-FUNCTIONAL THEORY FOR CRYSTALSPOSSESSING INTERNAL DEGREES OF FREEDOM戸田 直大 1) ，君塚 肇 2) ，尾方 成信 3)
Naohiro TODA, Hajime KIMIZUKA and Shigenobu OGATA
1) 大阪大学大学院基礎工学研究科
住友電気工業（株）(〒 560-8531

豊中市待兼山町 1-3,

E-mail: n-toda@mbox.me.es.osaka-u.ac.jp)(〒 664-0016

伊丹市昆陽北 1-1-1,

E-mail: toda-naohiro@sei.co.jp)

2) 大阪大学大学院基礎工学研究科(〒 560-8531

豊中市待兼山町 1-3,

E-mail: kimizuka@me.es.osaka-u.ac.jp)

3) 大阪大学大学院基礎工学研究科(〒 560-8531

豊中市待兼山町 1-3,

E-mail: ogata@me.es.osaka-u.ac.jp)

We proposed a new ﬁnite-element analysis method based on a numerical stress-strain
database to calculate the non-linear large elastic deformation for single crystals possessing internal degrees of freedom. The database is constructed during ﬁnite-element analysis
on-the-ﬂy based on the ﬁrst-principles density functional theory, which is able to obtainhigh-precision stress without any empirical parameters even under ﬁnite strained conditions. The database signiﬁcantly improves the total computational eﬃciency without
loss of accuracy. We also carried out parallel-computational method for the database
construction process to realize further improvement of computational eﬃciency. The effectiveness of the proposed inverse analysis method is demonstrated through numericalsimulation for several example problems.
Key Words : Quasicontinuum Method, Stress-Strain Database, Finite Element Method,
Parallel Computation, Diamond Crystalや応力状態によって大きく変化する (2) ．ナノレベルの実験

1. 緒言
ナノテクノロジーに注目が集まる中，ナノインデンテーションに代表されるナノスケールでの力学特性を測る実験技を理解するための計算機シミュレーションではこれらをどこまで正確に評価できるかが重要である．

術の向上はめざましい．それに伴って，これらのナノスケー

このようなナノレベルの解析対象に対して最も有効な力

ルの力学実験の詳細を理解するために，計算機シミュレー

学解析手法は，原子を直接扱う分子動力学法である．ただションが用いられるようになり，計算機シミュレーションにし，経験的な原子間相互作用を用いる古典分子動力学法でもナノスケールでの高精度な力学解析が求められている．多は，10 %もの変形に対して，経験的な原子間相互作用が有効結晶材料であっても，ナノインデンテーションのようにナノであるかどうかの疑問が常に残る．このような原子間相互作スケールの局所的な性質を調べる試験で結晶粒界以外を調べ用の不確かさのない方法として，密度汎関数法に代表されるた場合には，無欠陥の単結晶を測定していることとほぼ等価

非経験的な第一原理計算法がある．非経験的な第一原理計算となる．無欠陥の単結晶は，永久変形（塑性変形）をするま

法では，原子間の相互作用を，電子状態を解くことにより決でに，10 %以上も弾性的に変形することがある(1)(2)(3)．こ定するため，様々な材料の様々な変形状態に対して信頼できのレベルのひずみでは，材料はもはや線形弾性体として振舞る結果を得ることができる．また，非経験的に相互作用を計わず，弾性域であっても強い非線形性を呈することになる．

算できることから，実験等から得られる知見を用いることなまた，弾性域を超えて永久変形へと移行する力学条件も材料く計算が可能であり，予測的な計算結果が得られることも大きな特徴である．しかしながら，第一原理計算法は，このよ2007 年 12 月 29 日受付，2008 年 1 月 26 日受理

うな利点と引き換えに，計算量が古典分子動力学法に比べて極めて多いという欠点を持っている．その結果，現実的に扱える原子数が少なくなり，モデルの大きさや境界条件に大きな制約がある．
一方，構造解析の手段として広く用いられている有限要素は，データベースを構築する際に並列処理を採用する．

2. 第一原理密度汎関数法に基づく局所擬似連続体法
連続体では，物体点の変形は，次式によって各物体点で定義される変形勾配F によって一義的に表現できる．

F (x) =法は，本来多数の原子で構成されている材料を連続体として近似することにより，解析する際の自由度を原子団の自由度3N （N ：原子数）から，要素を結びつける節点の自由度 3M（M ：節点数，M ≪ N ）に減らし，大幅な高速化を実現している手法とみなすことができる．そして，境界条件の自由度も高い．しかし，ナノスケールの非線形弾性問題を解くにあたり，有限要素法には 2 つの大きな問題点がある．一つは，一般的な連続体を対象とした計算ではスケールの概念がなく，原子スケールの現象の詳細を獲得することはできない点である．ただし，単結晶中の数ナノメートルを超える長さスケールをもつ現象では，連続体近似が有効であるとの解析もあり (4) ，ナノインデンテーションを始めとするナノスケー
ルの計算の多くに適用されている (5) ．もう一つは，幾何学∂x∂X(1)

ここで，xとX はそれぞれ現配置と基準配置での物体点の位置ベクトルを表す．外部からの力学的負荷に対する材料の力学応答は，力学的負荷によって引き起こされた変形に起因する材料内部の自由エネルギーの変化によって記述できる．固体結晶材料の温度一定の試験ではエントロピーの寄与が無視できるので，自由エネルギーの変化は主に材料を構成する原子間のポテンシャルエネルギーの変化に起因している．転位やき裂等の欠陥の生成・消滅のない弾性変形をする系では連続体近似を適用することができ，材料全体のひずみエネル
ギーは次式に示すように，変形勾配分布の汎関数として表現できる．

的な変形状態と力学的な応力状態を結びつける構成式が必要

strainEtot≃ ∆U (F (x))

な点である．これまで材料の構成式は，その材料に対する力

ここで U は系に含まれる原子間のポテンシャルエネルギー学実験から得られる知見に基づいて，経験的に決定されていであり，ここで，変形勾配が零の状態を結晶の平衡状態であたのが現状である．この方法はマクロな解析では極めて有効るとし，そこを基準にポテンシャルエネルギーの変化を見積であるが，ナノスケールの構造解析に適用しようとすると，もっている．一方，材料全体のひずみエネルギーは，各物体10 %を超えるひずみまであらゆる変形状態を実現したナノ点 x0 で計算されるひずみエネルギー密度 e(x0 ) の積分で表スケール実験を実施し，そこから構成式を構築する必要があ現できる．

る．しかし，実際このような実験は，現状の実験技術をもってしても極めて実行困難である．特に，非線形弾性の振る舞いをすべて考慮するのは不可能である．

(2)

∫ ∫

strainEtot=VΩ(r ̸=0)e(F (x0 ), F (x), r)dx0 dr(3)

ここで V は材料の体積を表す．また，r は x0 から x へのベ

そこで，我々は，これら第一原理計算法と有限要素法それクトルである．e(x0 ) は一般的にその物体点 x0 での変形勾ぞれの長所を活かし，短所を互いに補う手法を開発した．構配だけでなく，物体点周りの領域 Ω 中の変形勾配にも依存成関係（応力–ひずみ関係）を非経験的で高精度な第一原理する．Ω の大きさは，均一な結晶の場合，原子間相互作用の密度汎関数計算手法に基づいて求めながら，有限要素法を実有効範囲によって決まる．尾方らの解析では，この有効範囲行する擬似連続体有限要素法である (4) ．

は面心立方金属で 2 ∼ 3 原子層，体心立方金属で 8 ∼ 9 原子

Fago ら (6) や，Hayes ら (7) は，計算量の少ない Orbital層程度と見積もられている (4) ．本解析で対象とするダイヤFree 密度汎関数法を用いた，同様の手法を提案している．しモンドは単一元素からなるため長距離のイオン性結合を持かしながら，Orbital Free 密度汎関数法は，電子密度の運動たず，かつ電子が局在した共有結合であるため，有効範囲はエネルギー汎関数の正確な形が不明であり，通常の密度汎関

高々2 ∼ 3 原子層である．例えば，ダイヤモンドの解析に使数に比べて精度が劣る．これらの簡易な密度汎関数法が利用用されている，経験的 Tersoﬀ ポテンシャル (8) では，ポテンされている背景は，通常の密度汎関数法では計算時間が長シャルのカットオフ距離を 3 Å 程度にしていることからも理く，有限要素計算全体の計算時間が現実的でなくなるためで解できる．よって，3 Å に対して十分大きなスケールの現象ある．そこで我々は，通常の密度汎関数法の精度を保ちつつ，を解析対象とする場合には，物体点のひずみエネルギーは，現実的な計算量の有限要素計算を実現するために，応力–ひその物体点の変形勾配のみの関数と考えても差し支えない．

ずみ関係の学習型データベースを用いる方法を提案し，これこれをここでは局所連続体近似と呼ぶ．
∫e(F (x0 ), F (x), r)dr ≃ e(F (x0 ))Ω(r ̸=0)を用いてアルミニウム単結晶の変形解析を実施した (4) ．
アルミニウムは基本格子に原子をひとつだけしか含まな(4)いため内部自由度を持たないが，本論文ではこの手法をダイ

この近似の下では次式のように，各物体点での変形勾配に起ヤモンドのように基本格子内に 2 個以上の原子を含む内部自因する物体点のポテンシャルエネルギーの増加分 ∆U (F (x0 ))由度のある系に適用できるよう拡張する．内部自由度のあるによって全体のひずみエネルギーが記述できる．
∫∫
strain∆U (F (x0 ))dx0e(F (x0 ))dx0 =Etot≃系の計算では計算量が内部自由度のない系に比べて多くなる
ため，計算のさらなる高速化が必要である．そこで本手法でVV(5)用いる (15) ．ABAQUS には，ユーザサブルーチンによって，z [001]構成関係式をカスタマイズできる機能がある．本解析ではこのユーザサブルーチンを作成し，本手法を ABAQUS 上で実現する．作成するユーザサブルーチンでは，入力される変形勾配F に対して，応力ﬀ と接線係数 J = ∆ﬀ/∆› を出力する．
本手法では，与えられた変形勾配に対する応力のデータを，aP3

aP1

後に説明する第一原理計算法を用いて構築されるデータベースから検索し，さらには接線係数を算出した上で，これらを

y [010]aP2合わせてサブルーチンの出力とする．

Primitive cell

まず，本サブルーチンには，有限要素プログラムから，要素サブルーチン内での具体的な計算は以下のとおりである．

x [100]Fig. 1

Crystal coordinates and primitive cell of diamondの変形に対応する変形勾配F が入力される．サブルーチン内で，変形勾配F を，回転テンソルRと右ストレッチテンソル

lattice.

U に極分解する．
F = RU

本手法では Cauchy-Born 則 (9) を適用し，物体点の変形を(7)原子配列のアフィン変形に置き換えることによってポテンシャ第一原理計算を行う座標系を固定し，モデルの結晶方位を任ルエネルギーの増加分を求める．Cauchy-Born 則の適用は，結意に選べるようにするためには，モデルの座標系から第一原

P′
晶の基本単位格子ベクトル aPn (n = 1, 2, 3) を an (n = 1, 2, 3)

理計算を行う座標系へ回転Qを施す必要がある．

へ変形勾配を用いて次式で変換することに対応する．このよU ′ = QU QT(8)

うに局所連続体近似の下で連続体と結晶構造を対応づける方
この変換後の右ストレッチテンソル U ′ に対応する応力 ﬀ ′ (U ′ )法をここでは局所擬似連続体法と呼ぶ．
PaP′
n = F anを，第一原理計算法によって構築した応力のデータベース(6)

か ら 求 め る ．そ し て ，得 ら れ た ﬀ ′ (U ′ ) に 対 し ，座 標 変 換

本解析では，この結晶のアフィン変形によるポテンシャルエ

ﬀ = RT QT ﬀ ′ (U ′ )QR を施し，モデルの座標系上での応力ネルギーの増分を第一原理密度汎関数法によって見積もる．と接線係数に変換した上で，有限要素プログラムに引き渡す．

また，局所擬似連続体に対する数値力学解析は有限要素法に3.3. 学習型データベースの構築よって実現する．

本手法で最も計算時間を要するのが各ストレッチテンソルに対する応力の第一原理計算である．そこで，計算の高速化

3. 計算手法のために，一度求めた応力はデータベースとして蓄積し，繰

3.1. 第一原理密度汎関数法計算

り返し利用することを考える．2 節に記したように，対称テ第一原理計算法として，一般化密度勾配近似 (10)（GGA：ンソルである右ストレッチテンソルU の 6 つの独立な成分をPerdew-Wang らにより提唱された形）された密度汎関数法

離散化し，6 次元の離散空間メッシュを作成する．そして，そを用いる．電子の波動関数を平面波により展開し，電子と原のメッシュ点それぞれに第一原理計算によって評価された応子の相互作用はウルトラソフト擬ポテンシャル (11) を用い力 6 成分を割り当て，これをデータベースとする．実際に必て表現する．具体的には計算コード VASP (Vienna Ab-initio

要な応力は離散化した点からの線形補間によって求める．本Simulation Package (12)(13) ) を用いて計算を行う．以下に示研究では，メッシュ間隔として，1 %を採用する．線形補間

すダイヤモンドに対する計算では，ブリルアンゾーン内の積には，各次元ごとに隣接するメッシュ上の応力が必要となる分を 5×5×5 の Monkhorst-Pack のk点メッシュ (14) によってため，26 個 (= 64 個) の隣接したメッシュ点上の応力を計算実行し，平面波のカットオフエネルギーを 358.2 eV と設定しておく必要がある．以前の計算で既に利用されたメッシュする．また，計算対象モデルであるスーパーセルには炭素原点では応力が既にもとまっており，これがそのまま利用でき子 2 個を含むダイヤモンドの基本単位格子を採用する（図 1るため，第一原理計算はこの 64 個のメッシュ点の内，いまだ参照）．本計算条件で，計算される格子定数の実験値との誤計算に利用されていないものに限り新たに実行することで足差は 1 %程度，せん断弾性定数も 3 %程度と極めて高い精度りる．

を実現している．本第一原理計算法を用いて，有限要素法で解析の初期段階では，データベースの蓄積が不十分であ与えられる変形勾配から単位体積あたりのポテンシャルエネり，多くの第一原理計算を必要とするが，解析が進むにつれルギー増加を求め，さらにそれから内部応力を計算する．応

データベースが充実し，必要とする第一原理計算の回数が減力を定義する座標系は図 1 に示した結晶座標系にとる．

り，その結果，計算時間が短縮される．また，データベース

3.2. 有限要素計算が充実した状態では，非線形有限要素計算で必要な繰り返し有限要素解析には，ABAQUS/Standard (V6.4.1, 2004) を収束計算や，有限要素計算での要素数が増えても，同じデータを何度も参照することになるため，必要とする第一原理計

search算の回数はほとんど増加しない．

MasterDatabase

データベース化のもうひとつの利点は，それを主としてrequest構築した研究者だけでなく，同時に複数の研究者がアクセスできる場所に置くことで，グループとしての作業効率を上げ

Job Scheduler

ることが可能な点にある．例えば，得られたデータベースをVASP jobs (0-512)

集中管理するサーバに置き，複数の研究者によって構成され...

るグループのメンバー全員が，参照およびデータ登録できる

Subtask Pool

環境を作る．これにより，データベースへの登録の計算を分散することができるとともに，他の研究者が計算し登録したデータを，別の研究者がリアルタイムに利用することが可能

reply

requestとなる．

3.4. 学習型データベースのデータ量の削減

Slave 1

Slave 2多くの結晶は，対称性を有するため，これに基づくデータ... Slave nベースの削減が可能である．ダイヤモンドのように，基本格子が立方構造の場合，3 つの基本格子ベクトルは等価であり，x，y ，z の三軸を交換してできる応力–ストレッチテンソル関

Fig. 2

Parallel task execution based on the master-slavemodel.

係はすべて等価である．この関係を利用し，以下に示すダイ
ヤモンドに関する計算では必要なデータベースを 1/6 に削減

い高速化を図ることが可能である．このようなタスク単位にしている．

よる粗粒度の並列化においては，通信量が相対的に小さく抑3.5. 内部自由度がある系への拡張えられるため並列性能にとっては有利に働く一方，負荷均衡

例えばダイヤモンドは，面心立方構造を持つ 2 つの副格（ロードバランス）に偏りが見られる場合はプロセッサ全体子によって構成されており，その相対位置ベクトルは，基本が有効に活用されずに並列化効率が抑制されることが多い．

単位格子を規定する 3 つのベクトルを基底ベクトルとする

特に，第一原理計算による構造最適化においては，最安定状と (0.25, 0.25, 0.25) であり，格子全体のひずみとは独立で態に収束するまでに要する反復数は初期値に強く依存するこある．よって，あるひずみを与えた状態で，最も自由エネルとから，各ジョブの負荷にはばらつきが存在する．このことギーが低くなるように，この相対位置（内部自由度）を最適は，負荷の少ない処理を担当したプロセッサに待ち時間が生化する必要がある．その際問題なのが，最適化計算の初期値じ，全体として並列性能が低減する要因となる．

依存性である．最適化計算を始める初期相対位置ベクトルにそこで本研究においては，ジョブスケジューラを活用した

よっては，最安定状態でない準安定状態に収束してしまう可タスクプール方式の並列化を採用することで，負荷分散を動能性がある．これをできるだけ回避するため，初期相対位置的に行い，各プロセッサにおける負荷均衡を保つことに努め

ベクトルを，基本セル内でいくつか用意し，そこから自由エた（図 2）．この方式では，必要となる 64 点の応力 ﬀ の内，ネルギー最適化計算を実行し，最も低い自由エネルギーに収データベースに登録されていない ﬀ のみを対象として，解析束した場合の応力値を，データベースに登録する．ダイヤモ

管理プロセッサ（マスター）が第一原理計算ジョブ群を単一ンドの場合では，この初期相対位置ベクトルに以下に示す 8キューのタスクプールに挿入する．そしてジョブスケジューケースを考えた．(0.25, 0.25, 0.25)，(0.29, 0.29, 0.29)，(0.21,ラが空いたプロセッサに対してジョブを動的に割り当てるこ

0.21, 0.21)，(0.29, 0.29, 0.13)，(0.46, 0.46, 0.46)，(0.42, 0.50,とにより，全体としてマスター・スレーブ型の並列処理が実0.42)，(0.54, 0.46, 0.38)，(0.46, 0.38, 0.46)．ただし，基底ベ現される．本方式では，各ジョブの負荷およびプロセッサ性

クトルは基本単位格子を規定する 3 つのベクトルである．な能が不均一であっても負荷分散が自動的に行われる利点があ

お，この 8 ケースの相対位置ベクトルは，事前の様々なひずり，間断なく，空いたプロセッサで必要な第一原理計算を実み下の第一原理原子構造安定化計算結果から経験的に決定行させることが可能となる．上述のように，結晶が内部自由した．

度を持つ場合は，1 組の応力–ひずみ関係を求めるに際して 83.6. 並列分散計算による高速化

種類の異なる初期値からの計算を行うため，最大 512 ケース

本手法での計算速度の律速因子は，データベース構築時に（= 64×8）の第一原理計算の試行を必要とする．スレーブ数必要な多数の第一原理計算であり，手法全体の高速化には，n が 512 に対して十分小さい場合では，各プロセッサに割りこの部分の高速化が最も効果的である．特に内部自由度を当てられる負荷のばらつきが平準化されるため，この並列分有する結晶の場合は，3.5 節に記したようにより多くの計算散処理手法を用いることで n の増加に対してほぼ線形に計ケース数を必要とすることから，この傾向は一層強まる．

算時間が短縮される．

これらの一連の第一原理計算に際しては，各ケースが独立に扱えることから，複数プロセッサを用いて並列実行を行

200

Axis 2 [111]Linear elastic theory

γ 12DFT
Present method

σ12 [GPa]

16012080
4000.0_Axis 3 [110]Fig. 3_Axis 1 [112]Calculation element and boundary conditions for adiamond crystal.

Fig. 40.10.2
γ120.30.4

Stress-strain relation under shear deformation of adiamond crystal. Triangles and squares represent the dataobtained from the present method and density-functionaltheory (DFT) calculations, respectively. For comparison, thedata from the linear elastic theory are indicated as circles.

4. 解析
4.1. ダイヤモンド単結晶の理想強度を応力–ひずみデータベースから線形補間で数値的に求めて以上の手法を用い，ダイヤモンド単結晶のすべり面 (111)いることによるものである．よって，計算時間を問わなけれ面内の [112] 方向に，せん断ひずみを与えた場合の解析を行ば，データベースの離散値の間隔を細かくするか，高次の補い，ダイヤモンドのせん断理想強度の評価を行った．

間法などの高精度な補間方法を用いることで，誤差を減らす

結晶の [112]，[111]，[110] 方向に 3 つの軸を取った立方体モデルを，同じ形状の 3 次元 8 節点 1 次要素 1 個で表現すことは可能である．

4.2. データベース化による高速化の検討

る．図 3 に示すように，この立方体の (111) 面の [112] 方向

同じく，ダイヤモンドによるせん断変形を通じて，データに，工学せん断ひずみ γ12 を与える．ここですべての節点でベース化による高速化の確認を行った．今回の計算で用いたせん断方向以外の変位を拘束する．

システムは，Linux 環境下において，デュアルコアを搭載し工学せん断ひずみ γ12 に対するせん断応力 σ12 の関係をRRた Intel⃝Pentium⃝
D Processor 840（3.20GHz）を 4 個使

図 4 に示す．図には本手法を用いた結果に加え，直接第一原用している（計 8 PE (Processing Element)）．データベース

理計算でせん断変形解析を実施した結果と線型弾性論に基づ更新，ジョブスケジューリング等の解析管理および有限要素く構成式を用いた結果も示している．本解析がデータベース計算（ABAQUS）のために 1 PE を使用し，残り 7 PE をデーを介してリファレンスとしている第一原理計算法の結果と，タベース構築のための第一原理計算（VASP）に割り当てた．

本手法による結果が一致しており，データベースが正常に機また，VASP コードや ABAQUS のユーザサブルーチンにつ能していることを示している．
本手法および第一原理計算法による結果は，ひずみが 5 %をR
いては，Intel⃝FORTRAN Compiler 7.1 によりオブジェクトコードに変換している．

超えたあたりから強い非線形性を呈するようになり，このひデータベースの構築前と構築後の状態で，それぞれ 4.1 節ずみ領域ではもはや線形弾性論が有効でないことがわかる．と同じ方向に 10 %のせん断ひずみを与えたモデルに対して

ただし，材料はいまだ弾性的に変形しており，せん断負荷を本手法を適用し，所要計算時間を調べた．まず，データベー除くと，もとの形状を回復する．

スによる高速化のみを評価するために，第一原理計算に割りせん断変形をさらに続けると，応力はやがて最大値に達し，当てる PE 数を 1 としている（図 5）．データベース構築前

その後減少を始める．これは，永久変形が始まったことを意において，応力値を第一原理計算から都度取得する場合は味しており，その開始応力を理想最大せん断応力と定義する．

31,683 秒を要するのに対し，データベース構築後は同解析

本手法を用いた場合，理想最大せん断応力は，約 32 %のひが 42 秒に短縮され，754 倍の高速化が図られることが確認

ずみで得られており，その値は 110 GPa である．直接第一原できた．ここで両者の差は，応力–ひずみデータベースを構

理計算による結果 101 GPa と比較すると，約 9 %の誤差で築するために要する第一原理計算群の処理時間に相当する．

ある．このように理想強度を非経験的に評価できるのが本手

3.6 節で述べたように，1 組の応力–ひずみ関係を求めるに際法の大きな特徴のひとつである．当然ではあるが線形弾性論して 512 ケースの第一原理計算を必要とすることから，本過ではそれが評価できない．この誤差は，本手法で得た応力–程において第一原理計算 1 ケース当たりの所要時間は 61.8

ひずみ曲線の微小なばらつきに起因している．これは，応力秒と見積もられる．それに対して応力–ひずみデータベース

35000(3) M.I.Eremets, I.A.Trojan, P.Gwaze, J.Huth, R.Boehler,
31683

V.D.Blank: The strength of diamond, Appl. Phys. Lett.，

Calculation time [s]30000

87(2005)， pp. 141902-1–141902-3.
25000(4) 尾方成信，磨矢一晃，渋谷陽二: 第一原理密度汎関数法

20000に基づく局所擬似連続体法, 日本機械学会論文集 A 編，

15000

71(2005)， pp. 1653–1658.
(5) 例 え ば ，T.Zhu, J.Li, K.J.Van Vliet, S.Ogata, S.Yip,10000

and S.Suresh: Predictive Modeling of Nanoindentation5000

Induced Homogeneous Defect Nucleation in Copper, J.
420
Before constructing
database

Fig. 5After constructing
database

Calculation time for a shear deformation analysis,

before and after constructing the stress-strain database.

Mech. Phys. Solids， 52(2004)， pp. 691–724.
(6) M.Fago, R.L.Hayes, E.A.Carter, and M.Ortiz: DensityFunctional-Theory-Based

Local

Quasicontinuum

Method: Prediction of Dislocation Nucleation, Phys.
Rev. B， 70(2004)， pp. 100102-1–100102-4.
(7) R.L.Hayes, G.Ho, M.Ortiz and E.A.Carter: Predictionを利用する場合は，データベースが格納されている記憶装置of Dislocation Nucleation during Nanoindentation ofへのアクセスと補間計算に要する時間のみで済むため，第一Al3 Mg by the Orbital-Free Density Functional Theory

原理計算と比較して無視できるほどの所要時間となる．

Local Quasicontinuum Method, Phil. Mag.， 86(2006)，

本解析手法では，データベースへのヒット率が低く多数の第一原理計算の呼び出しを必要とする初期段階がオーバー

pp. 2343–2358.
(8) J.Tersoﬀ: Empirical Interatomic Potential for Silicon

ヘッドとなるが，これらの第一原理計算群については上述のwith Improved Elastic Properties,

並列分散処理手法によりタスク分割が行えるため，本オー

38(1998)， pp. 9902–9905.

Phys. Rev. B，

バーヘッドは PE 数を増加させることでほぼ線形に減少した．

(9) J.L.Ericksen: Phase Transformations and Material In-

データベースの構築の前後における速度比は 754 倍に上stabilities in Solids, (1984)，Academic Press, New York．

ることから，並列分散処理による効果以上に，学習を終えた(10) J.P.Perdew, J.A.Chevary, S.H.Vosko, K.A.Jackson,後のデータベースを利用することの利点は大きい．特に多数

M.R.Pederson,の要素からなる実際の有限要素計算モデルでは，隣接する要Molecules, Solids, and Surfaces:

素は類似のひずみ状態（応力状態）になっている場合が多く，

the Generalized Gradient Approximation for Ex-

ある要素のひずみ状態に対して構築されたデータベースが，change and Correlation,その他の要素の計算にも利用できる．そのため，非線形有限pp. 6671–6687.

要素解析における収束計算の初期段階から，データベース化

D.J.Singh,C.Fiolhais:Atoms,Applications ofPhys. Rev. B， 46(1992)，(11) D.Vanderbilt: Soft Self-Consistent Pseudopotentials inの効果が期待できる．

a Generalized Eigenvalue Formalism,

5. 結言54(1990)， pp. 7892–7895.

Phys. Rev. B，結晶内部自由度を有する単結晶の弾性非線形大変形問題(12) G.Kresse, J.Hafner: Ab-initio Molecular Dynamics forを高精度かつ現実的な時間で計算できる第一原理密度汎関数

Liquid Metals, Phys. Rev. B， 47(1993)， pp. 558–561.

法に基づく局所擬似連続体法を提案した．局所擬似連続体に(13) G.Kresse, J.Furthmüller: Eﬃcient Iterative Schemes

対する数値計算法としては有限要素法を適用し，高速化手段

for Ab-initio Total-Energy Calculations using a Plane-には学習型データベースの適用および並列化計算を取り入れ

Wave Basis Set, Phys. Rev. B， 54(1996)， pp. 11169–た．そして，この手法をダイヤモンド単結晶のせん断変形解11186.

析に適用し，解析精度ならびにその計算速度を評価した．

参考文献(14) H.J.Monkhorst, J.D.Pack: Special Points for BrillouinZone Integrations, Phys. Rev. B， 13(1976)， pp. 5188–5192.

(1) S.Ogata, J.Li, S.Yip: Ideal Pure Shear Strength of Aluminum and Copper, Science，298(2002)，pp. 807–811.
(2) S.Ogata, J.Li, N.Hirosaki, Y.Shibutani, S.Yip: Ideal
Shear Strain of Metals and Ceramics, Phys. Rev. B，70(2004)， pp. 104104-1–104104-7.

(15) ABAQUS Version 6.4, Hibbit, Karlson, and Sorensen
Inc, (2004).