JP2013253917A - 電気伝導特性予測方法及びプログラム - Google Patents

Abstract

【課題】 本発明の課題は、グラフェンシートを使った電子材料の電気伝導特性予測を高速に行うことを目的とする。
【解決手段】 上記課題は、コンピュータによって実行される電気伝導特性予測方法であって、該コンピュータが、バンド計算によって、前記記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力することにより達成される。
【選択図】 図３

Description

本発明は、グラフェンシートを使った電子材料の電気伝導特性を予測する技術に関する。

グラフェンは、例えば、黒鉛結晶中において炭素の六角形格子を構成するｓｐ^２結合をした炭素原子よりなる原子層であるが、散乱の効果を抑制できれば室温でも２０００００ｃｍ²Ｖ^-1ｓ^-1を超える非常に大きな電子移動度を達成可能であることから、グラフェンのシートを使って超高速電子装置を作製する研究がなされている。

グラフェンシートはその製造過程において原子空孔欠陥が導入される可能性がある。また、非特許文献1のように意図的に原子空孔欠陥を導入する場合もある。そのような欠陥を含むグラフェンシートの電気伝導特性を知ることは電子装置へ応用する上で非常に重要である。

電子材料の原子レベルでの電気伝導特性を調べる方法としては、非平衡グリーン関数法に基づいた量子伝導計算が従来用いられる。この方法は、（１）まず半無限のリード電極で特性を調べたい対象の系を挟み込んだ構造のモデルを用意する。（２）次に、印加電圧に対応するポテンシャル勾配をそのリード電極間に発生させ、セルフコンシステントにポテンシャル電子状態を収束させる。（３）それにより得られたグリーン関数を用い、ランダウア公式に従って、逆格子空間の各ｋ点における伝導方向の透過係数を求める。（４）これをすべてのｋ点について求めて合計し、電流値を計算する。印加電圧を変えて以上の（２）−（４）を繰り返すことで電流電圧特性（電気伝導特性）を得ることができる。

J. Bai et al., "Graphen nanomesh", Nature Nanotechnology, March 2010, vol.5, p.190

上述した非平衡グリーン関数法に基づいた量子伝導計算の方法では、特に（２）のセルフコンシステントにポテンシャルと電子状態を収束させるプロセス、及びグリーン関数を作るプロセスに高コストの計算量を要する。

例えば、原子空孔欠陥を含むグラフェンシートの電気伝導特性を調べるためには、少なくとも１００原子程度から成る構造モデルを用意する必要があるが、その場合に任意の印加電圧下での電流値を十分な精度で計算するためには、例えば約１００ＭＰＩ並列の計算機を使って約２時間を要する。

よって、本発明の目的は、様々な欠陥パターンによる電気伝導特性の違いを探索したい場合や、更に多くの原子数を必要とする構造モデルを用いたい場合には、上述したような計算方法より低コストで電気伝導特性を予測することである。

開示の技術は、コンピュータによって実行される電気伝導特性予測方法であって、該コンピュータが、バンド計算によって、前記記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する。

また、上記課題を解決するための手段として、コンピュータに電気伝導特性を予測させるためのプログラム、そのプログラムを記録した記録媒体、及び、電気伝導特性予測装置とすることもできる。

開示の技術では、グラフェンシート中に原子空孔欠陥を含む種々の構造に関して、グラフェンシートの電流電圧特性を関連技術より低コストの計算で電気伝導特性を予測することができる。

コンピュータ装置のハードウェア構成を示す図である。 電気伝導特性を予測するための機能構成例を示す図である。 電気伝導特性予測処理における動作例を示す図である。 手順（２）での処理を説明するためのフローチャート図である。 手順（３）及び手順（４）での処理を説明するためのフローチャート図である。 手順（５）及び手順（６）での処理を説明するためのフローチャート図である。 第１の実施例において調査したい欠陥を含むグラフェンシートを用いた電子装置の模式図である。 第１の実施例における欠陥が二次元周期で並んだ構造を示す平面図である。 グラフェンの伝導方向のバンド構造を説明するための図である。 透過係数の比較を示す図である。 電流電圧特性の計算結果の比較例を示す図である。 第１の実施例において、二次元周期構造の伝導方向周期ｎ_Ａに対する電流電圧特性の依存性を示す図である。 第１の実施例において、伝導方向に垂直な方向のｋ点数ｎｋ_Ｚに対する電流電圧特性の依存性を示す図である。 第２の実施例に係る、調査したい欠陥を含むグラフェンシートを用いた電子装置の模式図である。 第２の実施例における電流電圧特性の計算結果を示す図である。

以下、本発明の実施の形態を図面に基づいて説明する。本実施の形態に係る、原子空孔欠陥を含むグラフェンシートの電気伝導特性予測の方法について、以下に箇条書きで順を追って説明する。

（１）調べたい（解析対象の）原子空孔欠陥を含むグラフェンシートについて、二次元周期で欠陥が並んだ構造モデルを作成する。

（２）バンド計算を行い、逆格子空間中の任意のk点において伝導方向のバンド構造を求める。

（３）フェルミ準位Ｅ_Ｆ前後０．８ｅＶ、即ち、

程度のエネルギー範囲にある各バンドについて傾きｄＥ／ｄｋの最大値を求め、ｄ_ｂ（Ｅ）とする。

（４）バンド計算により予め求めておいた原子空孔欠陥のない理想的なグラフェンシートのフェルミ準位Ｅ_Ｆ前後０．８ｅＶ、即ち、

のバンドの傾きｄ_ａと手順（３）で求めたｄ_ｂ（Ｅ）を用いて、透過係数Ｔ_ｋ（Ｅ）をＴ_ｋ（Ｅ）＝ｄ_ｂ（Ｅ）／ｄ_ａの式に従い算出する。

ここで、発明者は、バンド計算によって得られる理想的なグラフェンシートのバンドがＥ_Ｆ付近で直線的であり傾きが一定であることと、関連技術の方法で計算した理想的なグラフェンシートの透過係数がＥ_Ｆ付近で１となることに着目している。

伝導方向に垂直な方向のｋ点全てについて（２）から（４）の手順を行う。

（５）伝導方向に垂直な方向のｋ点全てについて（２）から（４）の手順を行った後、得られた透過係数Ｔ_ｋ（Ｅ）を合計してｋ点数ｎｋ_Ｚで規格化する。

（６）Ｔ_ｔｏｔ（Ｅ）を任意の印加電圧に対応するエネルギー範囲で積分して電流値を計算し、電流電圧特性を得る。

上述した手順（２）から手順（４）により、任意の構造の欠陥を含むグラフェンシートの電流電圧特性を計算することが可能となる。手順（１）から手順（４）で透過係数を求める方法が新規技術である。求めた透過係数を用いて、手順（５）、手順（６）で電流電圧特性を計算すること、また、手順（１）、手順（２）を透過係数の計算に適用することは知られていない。

上記手順（１）から手順（６）による原子空孔欠陥を含むグラフェンシートの電気伝導特性予測処理は、図１に示されるようなハードウェア構成を有するコンピュータ装置１００によって実行される。

図１は、コンピュータ装置のハードウェア構成を示す図である。図１において、コンピュータ装置１００は、コンピュータによって制御される端末であって、ＣＰＵ（Central Processing Unit）１１と、主記憶装置１２と、補助記憶装置１３と、入力装置１４と、表示装置１５と、出力装置１６と、通信Ｉ／Ｆ（インターフェース）１７と、ドライブ１８とを有し、バスＢに接続される。

ＣＰＵ１１は、主記憶装置１２に格納されたプログラムに従ってコンピュータ装置１００を制御する。主記憶装置１２には、ＲＡＭ（Random Access Memory）等が用いられ、ＣＰＵ１１にて実行されるプログラム、ＣＰＵ１１での処理に必要なデータ、ＣＰＵ１１での処理にて得られたデータ等を格納する。また、主記憶装置１２の一部の領域が、ＣＰＵ１１での処理に利用されるワークエリアとして割り付けられている。

補助記憶装置１３には、ハードディスクドライブが用いられ、各種処理を実行するためのプログラム等のデータを格納する。補助記憶装置１３に格納されているプログラムの一部が主記憶装置１２にロードされ、ＣＰＵ１１に実行されることによって、各種処理が実現される。記憶部１３０は、主記憶装置１２及び／又は補助記憶装置１３を有する。

入力装置１４は、マウス、キーボード等を有し、ユーザがコンピュータ装置１００による処理に必要な各種情報を入力するために用いられる。表示装置１５は、ＣＰＵ１１の制御のもとに必要な各種情報を表示する。出力装置１６は、プリンタ等を有し、ユーザからの指示に応じて各種情報を出力するために用いられる。通信Ｉ／Ｆ１７は、例えばインターネット、ＬＡＮ（Local Area Network）等に接続し、外部装置との間の通信制御をするための装置である。

コンピュータ装置１００によって行われる処理を実現するプログラムは、例えば、ＣＤ−ＲＯＭ（Compact Disc Read-Only Memory）等の記憶媒体１９によってコンピュータ装置１００に提供される。即ち、プログラムが保存された記憶媒体１９がドライブ１８にセットされると、ドライブ１８が記憶媒体１９からプログラムを読み出し、その読み出されたプログラムがバスＢを介して補助記憶装置１３にインストールされる。そして、プログラムが起動されると、補助記憶装置１３にインストールされたプログラムに従ってＣＰＵ１１がその処理を開始する。尚、プログラムを格納する媒体としてＣＤ−ＲＯＭに限定するものではなく、コンピュータが読み取り可能な媒体であればよい。コンピュータ読取可能な記憶媒体として、ＣＤ−ＲＯＭの他に、ＤＶＤディスク、ＵＳＢメモリ等の可搬型記録媒体、フラッシュメモリ等の半導体メモリであっても良い。

また、コンピュータ装置１００によって行われる処理を実現するプログラムが、通信Ｉ／Ｆ１７を介して外部装置から提供されてもよい。或いは、外部装置へ該プログラムを提供し、後述される各処理は外部装置で実現されるように構成してもよい。通信Ｉ／Ｆ１７による通信は無線又は有線に限定されるものではない。

電気伝導特性を予測するための機能構成について、図２及び図３で説明する。図２は、電気伝導特性を予測するための機能構成例を示す図である。図３は、電気伝導特性予測処理における動作例を示す図である。

図２において、コンピュータ装置１００は、入力パラメータ取得部４１と、構造モデル作成部４２と、バンド構造算出部４３と、透過係数算出部４４と、電流電圧特性算出部４５とを有する。入力パラメータ取得部４１と、構造モデル作成部４２と、バンド構造算出部４３と、透過係数算出部４４と、電流電圧特性算出部４５とは、ＣＰＵ１１が対応するプログラムを実行することによる処理によって実現される。

また、コンピュータ装置１００の記憶部１３０には、入力パラメータ５１と、構造モデル５２と、規定値５３と、出力データ５４とが格納される。

入力パラメータ取得部４１は、電気伝導特性予測処理を行う際に必要となる入力パラメータ５１を、ユーザから取得する処理部である。入力パラメータ取得部４１は、入力パラメータ５１を取得するために、入力パラメータ５１を設定可能なユーザインタフェースを表示装置１５に表示する。ユーザから取得した入力パラメータ５１は、記憶部１３０に格納される。

構造モデル作成部４２は、記憶部１３０に格納された入力パラメータ５１に基づいて、構造モデル５２を作成する処理部である。入力パラメータ５１に基づいて、ユーザが調べたい（解析対象の）原子空孔欠陥を含むグラフェンシートについて、二次元周期で欠陥が並んだ構造モデル５２が作成され、記憶部１３０に格納される。構造モデル作成部４２での処理は、上記手順（１）に相当する。

バンド構造算出部４３は、構造モデル５２を用いて、逆格子空間中の任意のk点において伝導方向のバンド構造を算出する処理部である。バンド構造算出部４３は、上記手順（２）に相当する。バンド構造算出部４３によって、原子空孔欠陥のない理想的なグラフェンシートのフェルミ準位Ｅ_Ｆ前後０．８ｅＶ）のバンドの傾きｄ_ａが算出される。

透過係数算出部４４は、透過係数Ｔ_ｋ（Ｅ）を算出する処理部であり、上記手順（３）及び（４）に相当する。透過係数算出部４４は、フェルミ準位Ｅ_Ｆ前後０．８ｅＶ程度のエネルギー範囲にある各バンドについて傾きｄＥ／ｄｋの最大値ｄ_ｂ（Ｅ）を求める（手順（３））。そして、透過係数算出部４４は、バンド構造算出部４３によって算出されたバンドの傾きｄ_ａと、手順（３）で求めた傾きｄＥ／ｄｋの最大値ｄ_ｂ（Ｅ）とを用いて、透過係数Ｔ_ｋ（Ｅ）をＴ_ｋ（Ｅ）＝ｄ_ｂ（Ｅ）／ｄ_ａの式に従い算出する（手順（４））。

図３に示すように、伝導方向に垂直な方向のｋ点全てについて、ｋ点毎に、バンド構造算出部４３と透過係数算出部４４とによって（２）から（４）の手順が行われる。得られた透過係数Ｔ_ｋ（Ｅ）は合計され、ｋ点数ｎｋ_Ｚで規格化される。

電流電圧特性算出部４５は、規格化された透過係数Ｔ_ｔｏｔ（Ｅ）を、入力パラメータ５１で指定される印加電圧に対応するエネルギー範囲で積分して電流値を計算し、電流電圧特性を得る処理部である。電流電圧特性算出部４５は、上記手順（６）に相当する。電流電圧特性算出部４５による算出結果（電流密度Ｃｄｅｎ）を含む出力データ５４が記憶部１３０に出力され格納される。

入力パラメータ５１には、原子座標、アームチェア方向のユニットセルサイズｎ_Ａ（＝８０）、ジグザグ方向のユニットセルサイズｎ_Ｚ、伝導方向（アームチェア方向）のｋ点数ｎｋ_Ａ（＝１０）、伝導方向に垂直な方向（ジグザグ方向）のｋ点数ｎｋ_Ｚ（〜ｎ_Ｚ／３５）、電流電圧特性予測に用いるバンドの本数ｎｂ（〜ｎ_Ａ＊ｎ_Ｚ／１０）、調べたい印加電圧の刻みｄＶ_ｂ（＝０．１Ｖ）、及び、調べたい印加電圧の最大値Ｖ_ｂ＿ｍａｘ（＝１．５Ｖ）等が含まれる。ここで、（ ）内の値は推奨値であるが、ユーザによって変更可能である。

入力パラメータ５１のうち、原子座標、ユニットセルサイズｎ_Ａ及びｎ_Ｚ、ｋ点数ｎｋ_Ａ及びｎｋ_Ｚ等が、構造モデル作成部４２によって構造モデル５２が作成される際に、参照される。

入力パラメータ５１のうち、電流電圧特性予測に用いるバンドの本数ｎｂ、調べたい印加電圧の刻みｄＶ_ｂ、及び、調べたい印加電圧の最大値Ｖ_ｂ＿ｍａｘ等が、電流電圧特性算出部４５によって参照される。

また、規定値５３には、グラフェンのＣ−Ｃボンド長

、強束縛近似計算で用いるホッピング積分値「２．５５ｅＶ」、欠陥のないグラフェンのフェミル順位付近のバンドの傾きｄ_a＝１０ｅＶ Ｂｏｈｒ等が含まれる。規定値５３は、透過係数算出部４４による透過係数Ｔ_ｋ（Ｅ）を算出する際に参照される。

出力データ５４は、電流電圧特性算出部４５による算出結果を示し、電流密度Ｃｄｅｎ等の値を含む。

以下に、手順（２）から（６）に係る処理ついて、図４から図６で説明する。図４は、手順（２）での処理を説明するためのフローチャート図である。図５は、手順（３）及び手順（４）での処理を説明するためのフローチャート図である。図６は、手順（５）及び手順（６）での処理を説明するためのフローチャート図である。

図４及び図５において、伝導方向に垂直な方向（ジグザグ方向）の各ｋ点毎（ｉｋ_ｚ番目毎）に、ＣＰＵ１１は、バンド構造算出部４３を起動して、構造モデル５２におけるジグザグ方向のｋ点での手順（２）の処理を行い、更に、透過係数算出部４４を起動して、手順（３）及び手順（４）での処理を行う。図６において、ＣＰＵ１１は、電流電圧特性算出部４５を起動して、手順（５）及び手順（６）での処理を行う。

先ず、図４において、バンド構造算出部４３による処理が開始され、バンド構造算出部４３はｊにゼロを設定して初期化する（ステップＳ１１）。そして、バンド構造算出部４３は、強束縛近似計算によって、ｍ番目のバンドのｉｋ_Ａ番目の逆格子点（０≦ｋ≦πの範囲、全部でｎｋ_Ａ個）夫々について固有エネルギーｅｉｇを計算する（ステップＳ１２）。

バンド構造算出部４３は、計算した固有エネルギーｅｉｇがゼロより大であるか否かを判断する（ステップＳ１３）。固有エネルギーｅｉｇがゼロより大であると判断した場合、ｊを１インクリメントして（ステップＳ１４）、ｅ（ｊ，ｉｋ_Ａ）に固有エネルギーｅｉｇを設定する（ステップＳ１５）。そして、バンド構造算出部４３は、ステップＳ１６へと進む。

一方、ステップＳ１３において、固有エネルギーｅｉｇがゼロ以下であると判断した場合、バンド構造算出部４３は、最後のｍ番目のバンドであるか否かを判断する（ステップＳ１６）。最後のｍ番目のバンドでない場合、バンド構造算出部４３は、ｍを１インクリメントして、ステップＳ１２へと戻り、上記同様の処理を行う。

一方、ステップＳ１６において、最後のｍ番目のバンドであると判断した場合、バンド構造算出部４３による手順（２）での処理は終了し、手順（３）及び手順（４）による処理へと進む（図５）。

図５中、ステップＳ２３が手順（３）に相当し、ステップＳ２８が手順（４）に相当する。図５において、透過係数算出部４４による処理が開示され、透過係数算出部４４はバンド数ｉｂ及びｋ点数ｉｋ_Ａに１を設定する（ステップＳ２１）。

透過係数算出部４４は、傾きの最大値ｄ_ｂ＿ｍａｘ（ｉｂ）及びエネルギー値ｅｎｅ（ｉｋ_ｚ，ｉｂ）をゼロに初期化する（ステップＳ２２）。

そして、透過係数算出部４４は、現在のｋ点（ｉｋ_Ａ）と次のｋ点（ｉｋ_Ａ＋１）との距離ｄｋ_Ａと、傾きｄ_ｂとを下記式により計算する（ステップＳ２３）。

傾きｄ_ｂは、次のｋ点（ｉｋ_Ａ＋１）のエネルギー値ｅ（ｉｂ，ｉｋ_Ａ＋１）と現在のｋ点（ｉｋ_Ａ）のエネルギー値ｅ（ｉｂ，ｉｋ_Ａ）との差を距離ｄｋ_Ａで割ることにより求められる。

次に、透過係数算出部４４は、傾きｄ_ｂが最大値ｄ_ｂ＿ｍａｘ（ｉｂ）より大であるか否かを判断する（ステップＳ２４）。大である場合、透過係数算出部４４は、最大値ｄ_ｂ＿ｍａｘ（ｉｂ）に傾きｄ_ｂを設定し、ｅｎｅ（ｉｋ_ｚ，ｉｂ）にｅ（ｉｂ，ｉｋ_Ａ）を設定して（ステップＳ２５）、ステップＳ２６へと進む。

一方、ステップＳ２４において、傾きｄ_ｂが最大値ｄ_ｂ＿ｍａｘ（ｉｂ）以下であると判断した場合、透過係数算出部４４は、現在のｋ点（ｉｋ_Ａ）がｋ点数ｎｋ_Ａから１減算した値に等しいか否かを判断する（ステップＳ２６）。等しくない場合、透過係数算出部４４は、ｉｋ_Ａを１インクリメントして（ステップＳ２７）、ステップＳ２３へ戻り、上記同様の処理を繰り返す。

一方、等しい場合、透過係数算出部４４は、伝導方向に垂直な方向のｋ点（ｉｋ_Ｚ）における透過係数ｔｒを計算する（ステップＳ２８）。透過係数ｔｒは、最大値ｄ_ｂ＿ｍａｘ（ｉｂ）をグラフェンのフェミル順位付近のバンドの傾きｄ_aで割ることによって得られる。

そして、透過係数算出部４４は、バンド番号ｉｂがバンド総数ｎｂと等しいか否かを判断する（ステップＳ２９）。バンド番号ｉｂがバンド総数ｎｂ等しくない場合、バンド総数ｎｂに達していないため、透過係数算出部４４は、バンド番号ｉｂを１インクリメントして（ステップＳ３０）、ステップＳ２２へと戻り、上記同様の処理を繰り返す。

一方、バンド番号ｉｂがバンド総数ｎｂと等しい場合、伝導方向に垂直なジグザグ方向のｋ点番号ｉｋ_Ｚでの透過係数算出部４４による処理は、ＣＰＵ１１によって一旦終了する。ＣＰＵ１１は、伝導方向に垂直なジグザグ方向のｋ点番号ｉｋ_Ｚがジグザグ方向のｋ点の総数ｎｋ_Ｚ以下であるか否かを判断する（ステップＳ３１）。ジグザグ方向のｋ点番号ｉｋ_Ｚが総数ｎｋ_Ｚ以下である場合、ＣＰＵ１１は、ステップＳ１１（図４）へ戻り、上述したように、バンド構造算出部４３を起動して、次のジグザグ方向のｋ点での手順（２）の処理を行わせ、更に、透過係数算出部４４を起動して、手順（３）及び手順（４）での処理を行わせる。

一方、ジグザグ方向のｋ点番号ｉｋ_Ｚが総数ｎｋ_Ｚを超えた場合、ＣＰＵ１１は、電流電圧特性算出部４５を起動して、ステップＳ４１（図６）へと進む。

図６において、電流電圧特性算出部４５による処理が開始され、電流電圧特性算出部４５は、ｉＶ_ｂに１を設定し（ステップＳ４１）、電流値Ｃｕｒｒにゼロを設定して初期化する（ステップＳ４２）。

そして、電流電圧特性算出部４５は、エネルギー値ｅｎｅ（ｉｋ_ｚ，ｉｂ）がｄＶ_ｂ＊ｉＶ_ｂ／２以下であるか否かを判断する（ステップＳ４３）。ｄＶ_ｂは、記憶部１３０に格納されている入力パラメータ５１から取得する。ｉＶ_ｂは、印加電圧の刻み回数である。

エネルギー値ｅｎｅ（ｉｋ_ｚ，ｉｂ）がｄＶ_ｂ＊ｉＶ_ｂ／２より大きい場合、電流電圧特性算出部４５は、電流値Ｃｕｒｒに、ｉｋ_ｚでの透過係数とエネルギー値とに基づくバンド間のエネルギー値を付加して（ステップＳ４４）、ステップＳ４５へと進む。一方、ｄＶ_ｂ＊ｉＶ_ｂ／２以下の場合、電流電圧特性算出部４５は、ステップＳ４５へと進む。

電流電圧特性算出部４５は、バンド番号ｉｂがバンド総数から１減算した値と等しいか否かを判断する（ステップＳ４５）。等しくない場合、バンド番号ｉｂを１インクリメントして（ステップＳ４６）、ステップＳ４３へと戻り、上記同様の処理を繰り返す。

一方、等しい場合、電流電圧特性算出部４５は、ジグザグ方向のｋ点番号ｉｋ_ｚがジグザグ方向のｋ点総数ｉｋ_ｚに等しいか否かを判断する（ステップＳ４７）。等しくない場合、電流電圧特性算出部４５は、ジグザグ方向のｋ点番号ｉｋ_ｚを１インクリメントして（ステップＳ４８）、ステップＳ４３へと戻り、上記同様の処理を繰り返す。

一方、ジグザグ方向のｋ点総数ｉｋ_ｚに等しい場合、電流電圧特性算出部４５は、電流密度Ｃｄｅｎを算出する（ステップＳ４９）。

電流密度Ｃｄｅｎは、上記数１０によって求められる。

電流電圧特性算出部４５は、印加電圧Ｖ_ｂ（＝ｄＶ_ｂ＊ｉＶ_ｂ）を計算し、計算した印加電圧と、ステップＳ４９で算出した電流密度Ｃｄｅｎとを含む出力データを記憶部１３０に格納する（ステップＳ５０）。

そして、電流電圧特性算出部４５は、印加電圧Ｖ_ｂが、記憶部１３０に格納されている入力パラメータ５１の調べたい印加電圧の最大値Ｖ_ｂ＿ｍａｘより小さいか否かを判断する（ステップＳ５１）。最大値Ｖ_ｂ＿ｍａｘより小さい場合、電流電圧特性算出部４５は、印加電圧の刻み回数ｉＶ_ｂを１インクリメントして（ステップＳ５２）、ステップＳ４２へ戻り、上述同様の処理を繰り返す。一方、印加電圧Ｖ_ｂが最大値Ｖ_ｂ＿ｍａｘ以上を示す場合、ＣＰＵ１１は、電流電圧特性算出部４５による処理は終了し、この電気伝導特性予測処理全体を終了する。

上述した電気伝導特性予測処理による実施例を以下に説明する。

［第１の実施例］
図７は、第１の実施例において調査したい欠陥を含むグラフェンシートを用いた電子装置の模式図である。第１の実施例では、グラフェンシート１０Ｇ中に欠陥１０Ａが存在し（構造ｂ）、その欠陥１０Ａを図７の左右から挟み込む形で設けられたグラフェン電極１０Ｓとグラフェン電極１０Ｄとの間に任意の電圧が印加された場合の両電極間に流れる電流を、上述した手順（１）から（６）に従って計算する。

伝導方向は、図７の左右方向で、グラフェンシート１０Ｇのアームチェア方向である。なお、図７、及び、以降の図では、簡単のために欠陥１０Ａの周りのＣ原子は未結合手を終端せずに描いているが、第１の実施例中では未結合手は全て終端されているものとして計算している。

手順（１）において、構造モデル作成部４２によって作成された、二次元周期で欠陥１０Ａが並んだ構造モデルを図８に示す。手順（２）以降は、図８中に示したユニットセル１１について計算が行われる。後述の通り、このユニットセル１１の伝導方向、即ち、アームチェア方向の周期サイズｎ_Ａは精度を決定するパラメータであるが、ここではｎ_Ａ＝８０（長さ約３４ｎｍ）を用いる。図８のユニットセル１１はｎ_Ａ＝５である。

手順（２）のバンド計算の手法としては、例えば、上述したように強束縛近似法を用いる。第一原理法等のその他のバンド計算手法を用いてもよいが、計算コスト面で強束縛近似法を用いるのが最も有利である。強束縛近似計算で用いる基底は炭素原子当たり１つのｐ_ｚ軌道で、相互作用は最近接原子間のみを考慮する。第一原理計算から求めた理想的なグラフェンのバンド構造を再現するように、ホッピング積分の値は２．５５ｅＶを採用する。

図９は、グラフェンの伝導方向のバンド構造を説明するための図である。ユニットセル１１と同じセルサイズで欠陥１０Ａのない理想的なグラフェン（構造ａ）と、欠陥１０Ａを含むグラフェン（構造ｂ）のユニットセル１１とについて求めた、Ｋ（０，２／３π，０）点からＫＸ（π，２／３π，０）点までの伝導方向のバンド構造をそれぞれ図９（Ａ）及び図９（Ｂ）に示す。ここで用いている強束縛近似計算によって得られるグラフェンのバンド構造はＥ_Ｆを軸に対称である。

図９（Ａ）に示したように、理想的なグラフェン（構造ａ）は０＜Ｅ−Ｅ_Ｆ＜０．８のエネルギー範囲では傾きｄＥ／ｄｋが一定であり、その値をｄ_ａとする。第１の実施例の場合、ｄ_ａ＝１０ｅＶ Ｂｏｈｒである。

一方、図９（Ｂ）に示した、欠陥１０Ａを含むグラフェンは、Ｅ_Ｆ付近で傾きがゼロに近いフラットなバンドを持ち、エネルギーが上がるに連れて傾きｄ_ａに近いバンドを持つようになり、傾きの最大値ｄ_ｂを示すようになる。

図９（Ｂ）で得られた各バンドにおいて、傾きが最大となるエネルギー値Ｅ−Ｅ_Ｆとその時の傾きｄ_ｂ（Ｅ−Ｅ_Ｆ）を求める。これが手順（３）である。

次に、手順（４）として、透過係数Ｔ_ｋ（Ｅ−Ｅ_Ｆ）をＴ_ｋ（Ｅ−Ｅ_Ｆ）＝ｄ_ｂ（Ｅ−Ｅ_Ｆ）／ｄ_ａの式に従って算出する。

図１０は、透過係数の比較を示す図である。図１０中、縦軸に透過係数Ｔ_ｋ（Ｅ−Ｅ_Ｆ）を示し、横軸にエネルギーＥ−Ｅ_Ｆ（ｅＶ）を示す。図１０において、理想的なグラフェン（（構造ａ）Ｉｄｅａｌ）、本第１の実施例（（構造ｂ）Ｎｅｗ）、及び、本第１の実施例と同じ構造について関連技術を用いた場合（（構造ｂ）Ｒｅｌａｔｅｄ Ａｒｔ）の夫々の透過係数Ｔ_ｋ（Ｅ−Ｅ_Ｆ）の計算結果が示されている。理想的なグラフェン（（構造ａ）Ｉｄｅａｌ）は、図９（Ａ）で示されるバンドの傾きｄａを示している。

図１０から、本実施の形態に係る計算によって得られる透過係数は、関連技術である非平衡グリーン関数法に基づいた量子伝導計算による結果をよく再現することが分かる。第１の実施例において、図１０にＩｄｅａｌとして示したように理想的なグラフェン（構造Ａ）が０＜Ｅ−Ｅ_Ｆ＜０．８のエネルギー範囲でＴ_ｋ（Ｅ−Ｅ_Ｆ）＝１であることを利用している。

（構造ｂ）関連技術では、図１０の透過係数の算出結果を得るのにＡＭＤ Ｏｐｔｅｒｏｎ ２．４ＧＨｚ ４コア×２の構成のコンピュータ装置を使って１０４ＭＰＩ並列で１２０分を要しているが、第１の実施例では、ＡＭＤ Ｏｐｔｅｒｏｎ ３．０ＧＨｚ ２コア×２の構成のマシン（以降マシンＡと呼ぶことにする）を使って１ＭＰＩ並列で８分しかかかっていない。単純に見積もって、計算コストを（３．０＊１＊８）／（２．４＊１０４＊１２０）＝１／１２４８にできたと言える。

最後に、手順（５）及び手順（６）において、以上の手順から得た透過係数を使って、電流電圧特性が算出される。手順（５）で必要となる伝導方向に垂直な方向のｋ点数ｎｋ_Ｚは後述の通り精度を決定するパラメータであるが、ここではｎｋ_Ｚ＝７を用いる。手順（６）で、例えば印加電圧Ｖ_ｂ＝１Ｖのときの電流密度を求めたい場合は、

を

のエネルギー範囲で積分し、その積分値をｕｐスピンとｄｏｗｎスピンの両方を考慮して２倍する。

第１の実施例ではＥ_Ｆを軸に対称な透過係数であるから、０＜Ｅ−Ｅ_Ｆ＜０．５でＴ_ｔｏｔ（Ｅ−Ｅ_Ｆ）の積分を行い、その積分値を４倍した値をＶ_ｂ＝１Ｖのときの電流密度とする。

このようにして得られた電流電圧特性を図１１に示す。図１１は、電流電圧特性の計算結果の比較例を示す図である。図１１中、縦軸に電流密度（μＡ／ｎｍ）を示し、横軸にバイアス電圧Ｖ_ｂ（Ｖ）を示す。

図１１から、０＜Ｖ_ｂ＜１Ｖの範囲においては、本第１の実施例（（構造ｂ）Ｎｅｗ）は関連技術（（構造ｂ）Ｒｅｌａｔｅｄ Ａｒｔ）と誤差最大１割程度の精度で電流電圧特性を予測できていることが分かる。透過係数が求まってから後のプロセスに要する計算コストは、本第１の実施例と関連技術でほとんど差がない。すなわち、本第１の実施例は関連技術の約１／１０００の計算コストで図１１の結果を得ることができる。

以上のように、第１の実施例を用いることで、欠陥を含むグラフェンシートの電流電圧特性を予測する際に、計算精度をほとんど損なうことなく、大幅な計算コストの削減が可能である。

ここで、第１の実施例において精度を決定するパラメータｎ_Ａ及びｎｋ_Ｚについて説明する。

先ず、手順（１）で用意する図８中のユニットセル１１の伝導方向、すなわちアームチェア方向の周期サイズｎ_Ａを上から順に１６０、８０、４０と変化させたときの電流電圧特性を図１２に比較した。図１２は、第１の実施例において、二次元周期構造の伝導方向周期ｎ_Ａに対する電流電圧特性の依存性を示す図である。図１２中、縦軸に電流密度（μＡ／ｎｍ）を示し、横軸にバイアス電圧Ｖ_ｂ（Ｖ）を示す。

図１２より、ｎ_Ａが大きいほど、図１０におけるエネルギーの刻みを細かくすることができ、精度が良くなる。ただし、本方法では計算コストのほとんどはバンド計算において固有値の対角化に要する時間であるため、ｎ_Ａを大きくすると、（ｎ_Ａ）^３でスケーリングして計算コストが高くなる。

例えば、マシンＡで１ＭＰＩ並列の場合、ｎ_Ａ＝１６０、８０、４０に対して、ｋ点１つあたりの透過係数を求める手順（４）終了までの計算時間はそれぞれ６０分、８分、１分である。ｎ_Ａが２倍になると、約２^３＝８倍の計算コストになっていることが分かる。第１の実施例では精度と計算コストのバランスからｎ_Ａ＝８０を選択した。

次に、手順（５）で必要となる伝導方向に垂直な方向のｋ点数ｎｋ_Ｚを上から順に１３、７、５、３と変化させたときの電流電圧特性の比較を図１３に示す。図１３は、第１の実施例において、伝導方向に垂直な方向のｋ点数ｎｋ_Ｚに対する電流電圧特性の依存性を示す図である。図１３中、縦軸に電流密度（μＡ／ｎｍ）を示し、横軸にバイアス電圧Ｖ_ｂ（Ｖ）を示す。

図１３より、ｎｋ_Ｚが大きいほど精度が高いことが分かる。計算コストのスケーリングはｎｋ_Ｚに等倍である。透過係数を求める手順（４）までの計算はｋ点毎に独立して行えるため、簡単に並列処理できる。第１の実施例では精度と計算コストのバランスからｎｋ_Ｚ＝７を選択した。この場合、伝導方向の逆格子ベクトル上にｋ点が約０．０４Ｂｏｈｒ^−１の間隔で並ぶことになる。

［第２の実施例］
図１４は、第２の実施例に係る、調査したい欠陥を含むグラフェンシートを用いた電子装置の模式図である。図１４（Ａ）に構造ｃのグラフェンシート２０Ｇと、図１４（Ｂ）に構造ｄのグラフェンシート３０Ｇを示す。

図１４（Ａ）に示すグラフェンシート２０Ｇは、欠陥２０Ａが存在する構造ｃを有し、欠陥２０Ａを左右から挟み込むようにグラフェン電極２０Ｓ及び２０Ｄが設けられる。図１４（Ｂ）においても同様に、グラフェンシート３０Ｇは欠陥３０Ａが存在する構造ｄを有し、欠陥３０Ａを左右から挟み込むようにグラフェン電極３０Ｓ及び３０Ｄが設けられる。

第２の実施例では、図１４（Ａ）のグラフェンシート２０Ｇを挟むグラフェン電極２０Ｓと２０Ｄとの間、図１４（Ｂ）のグラフェンシート３０Ｇを挟むグラフェン電極３０Ｓと３０Ｄとの間に任意の電圧が印加された場合の両電極間に流れる電流を上述した手順（１）から手順（６）に従って計算する。

伝導方向は、図１４（Ａ）及び図１４（Ｂ）の左右方向で、グラフェンシート２０Ｇ及び３０Ｇのアームチェア方向である。なお、図１４（Ａ）及び図１４（Ｂ）では、簡単のために、欠陥２０Ａ及び欠陥３０Ａの周りのＣ原子は、未結合手を終端せずに描いているが、第２の実施例中では未結合手は全て終端されているものとして計算している。

精度を決定するパラメータとして、ｎ_Ａ＝８０、ｎ_ｋ＝７を用いて、第２の実施例と同様の方法で求めた電流電圧特性を図１５に示す。図１５は、第２の実施例における電流電圧特性の計算結果を示す図である。図１５中、縦軸に電流密度（μＡ／ｎｍ）を示し、横軸にバイアス電圧Ｖ_ｂ（Ｖ）を示す。

第１の実施例同様、第２の実施例での計算結果（Ｎｅｗ）は関連技術（Ｒｅｌａｔｅｄ Ａｒｔ）による計算結果をよく再現していることが分かる。計算コストも第１の実施例とほぼ同じで、手順（４）終了までの計算時間は例えばマシンＡの１ＭＰＩ並列でどちらも８分である。

上述したように、本実施の形態により、グラフェンシート中に原子空孔欠陥を含む構造について、関連技術による非平衡グリーン関数法に基づいた量子伝導計算を用いる方法と同程度の精度を保ちつつ計算コストを削減することができる。また、様々な欠陥パターンによる電気伝導特性の違いの探索や、数百原子以上から成る構造モデルを要する欠陥を含むグラフェンシートの電気伝導特性の調査も可能となる。

本発明は、具体的に開示された実施例に限定されるものではなく、特許請求の範囲から逸脱することなく、種々の変形や変更が可能である。

以上の第１から第２の実施例を含む実施形態に関し、更に以下の付記を開示する。
（付記１）
コンピュータによって実行される電気伝導特性予測方法であって、該コンピュータが、
バンド計算によって、前記記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、
所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、
原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、
前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、
所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する
ことを特徴とする電気伝導特性予測方法。
（付記２）
所定の原子空孔欠陥を有するグラフェンシートについて、二次元周期で欠陥が並んだ前記構造モデルを作成して、前記記憶部に記憶する
ことを特徴とする付記１記載の電気伝導特性予測方法。
（付記３）
前記構造モデルの伝導方向の周期サイズは３４ｎｍ以上である
ことを特徴とする付記２記載の電気伝導特性予測方法。
（付記４）
前記伝導方向に垂直な方向のｋ点数を、前記伝導方向の逆格子ベクトル上にｋ点が約０．０４Ｂｏｈｒ^−１の間隔で並ぶように定められている
ことを特徴とする付記１乃至３のいずれか一項記載の電気伝導特性予測方法。
（付記５）
バンド計算によって、記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、
所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、
原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、
前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、
所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する
処理をコンピュータに実行させるプログラム。
（付記６）
原子空孔欠陥を有するグラフェンシートの構造モデルを記憶した記憶部と、
バンド計算によって、前記記憶部に記憶された前記構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得するバンド構造算出部と、
所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出する透過係数算出部と、
前記伝導方向に垂直な方向のｋ点毎に、前記バンド構造算出部と透過係数算出部とによる処理を行って得られた透過係数を合計してｋ点数で規格化し、所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する電流電圧特性算出部と
を有することを特徴とする電気伝導特性予測装置。

１１ ＣＰＵ
１２ 主記憶装置
１３ 補助記憶装置
１４ 入力装置
１５ 表示装置
１６ 出力装置
１７ 通信Ｉ／Ｆ
１８ ドライブ
１９ 記憶媒体
４１ 入力パラメータ取得部
４２ 構造モデル作成部
４３ バンド構造算出部
４４ 透過係数算出部
４５ 電流電圧特性算出部
５１ 入力パラメータ
５２ 構造モデル
５３ 規定値
５４ 出力データ
１００ コンピュータ装置
１３０ 記憶部

Claims (5)

1. コンピュータによって実行される電気伝導特性予測方法であって、該コンピュータが、
   バンド計算によって、前記記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、
   所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、
   原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、
   前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、
   所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する
   ことを特徴とする電気伝導特性予測方法。
2. 所定の原子空孔欠陥を有するグラフェンシートについて、二次元周期で欠陥が並んだ前記構造モデルを作成して、前記記憶部に記憶する
   ことを特徴とする請求項１記載の電気伝導特性予測方法。
3. 前記構造モデルの伝導方向の周期サイズは３４ｎｍ以上である
   ことを特徴とする請求項２記載の電気伝導特性予測方法。
4. 前記伝導方向に垂直な方向のｋ点数を、前記伝導方向の逆格子ベクトル上にｋ点が約０．０４Ｂｏｈｒ^−１の間隔で並ぶように定められている
   ことを特徴とする請求項１乃至３のいずれか一項記載の電気伝導特性予測方法。
5. バンド計算によって、記憶部に記憶された原子空孔欠陥を有するグラフェンシートの構造モデルにおける逆格子空間中の任意のｋ点において伝導方向のバンド構造を取得し、
   所定のエネルギー範囲にある各バンドについて傾きの最大値を取得し、
   原子空孔欠陥のない理想的なグラフェンシートの前記所定のエネルギー範囲にあるバンドの傾きで該最大値を割ることによって透過係数を算出し、
   前記伝導方向に垂直な方向のｋ点毎に前記バンド構造と前記最大値とを取得して算出した前記透過係数を合計してｋ点数で規格化し、
   所定の印加電圧に対応するエネルギー範囲で積分することによって電流情報を取得して前記記憶部に出力する
   処理をコンピュータに実行させるプログラム。

Images