SENAC Vol. 48, No. 4（2015. 10）— 1 —［共同研究成果］

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション
菅野 学,1 小林倫仁,1 Wilfredo Credo Chung,1,2 山崎 馨,1,3 瀬高 渉,4 河野裕彦 1
1
東北大学大学院理学研究科化学専攻
2
Department of Chemistry, De La Salle University−Manila, Philippines
3
北海道大学大学院理学研究院化学部門
4
首都大学東京大学院都市環境科学研究科分子応用化学域

動的機能分子の 1 つとして有望視される分子ジャイロスコープは外部骨格（固定子）によって保護された環（回転子）を持ち, それが両者の間の化学結合を軸として回転する. 近年合成されたフッ素置換回転子を有する結晶性分子ジャイロスコープ（RotF）は, テラヘルツ波による内部回転の直接励起が可能と考えられている. その実現に向けて, 我々は密度汎関数法の 1/10 から 1/100の計算負荷で同程度の計算精度が期待できる密度汎関数強束縛（DFTB）法を用いて, RotF の結晶構造と回転動力学のシミュレーションを行った. まず, DFTB 法の計算プログラム DFTB+を並列コンピュータ NEC LX 406Re-2 に実装し, その並列化性能を検証した. 次に, 内部回転の有効ポテンシャルエネルギー曲線の探索と, 熱運動の分子動力学シミュレーションを実行した. 結果の詳細な解析から, RotF の活性化エネルギーの値や室温における回転頻度を予測した.

1. 序論
新規な動的機能を発現する分子の創製およびその制御はナノテクノロジーの中心的課題である.
ナノサイズの分子は単分子でも多数の自由度を有するが, それらがファンデルワールス力で結び付いた固体（分子性結晶）においては隣接する分子が互いに相関し, 超多次元系を形成する. その超多次元空間の中から如何にして有効な反応座標を抽出し, ポテンシャルエネルギー曲面の鞍点に相当する遷移状態を制御するかという問題は, 分子の動的機能を高効率で操作するための方法論と直結する重要な学術的テーマである.
21 世紀に入り, 米国の Garcia-Garibay らは外部骨格（固定子）によって保護された環（回転子）が結合を軸として回転する分子ジャイロスコープの概念を提案し, 実際に種々の結晶性分子ジャイロスコー
固定子プを合成した. 図 1 にその代表例と核磁気共鳴（Nuclear magnetic resonance; NMR）スペクトル測定
回転子から見積もられた回転の活性化エネルギーを示した
回転軸[1-3]. 嵩高い固定子ほど結晶中で回転子の周囲に大きな空間を確保し, 活性化エネルギーを低下（遷移状態を安定化）させる.
瀬高らは固定子の設計によって遷移状態の安定性を調節するこのアイディアを発展させ, 大きく頑強12.8 kcal/mol< 5 kcal/molなジシロキサンかご状骨格を持つ分子ジャイロスコープ（図 2 で X = H としたもの, 以下 RotH と表記） 図 1 Garcia-Garibay らが合成した結晶性分子ジャイロスコープの例. 数値はを合成した[4]. 結合軸の両端を架橋したジシロキサ実測された回転の活性化エネルギー
ン鎖が周囲の影響を軽減し, フェニレン環は室温で[1-3].
NMR の時間スケールよりはるかに速く回転する.

— 2 —

SENAC Vol. 48, No. 4（2015. 10）

河野らは RotH の結晶構造を理論計算で再現することに成功し, 活性化エネルギーが図 1 の従来型と比べて非常に低い（0.7 kcal/mol）ことを示した[5].
瀬高らが合成した分子ジャイロスコープの複屈折性が温度に依存して劇的に変化する様子が偏光顕微鏡観察により確認されている[6]. このことから, 環の内部回転を瞬時に誘起できれば, 液晶よりも高速に動作する光学材料としての応用が見込まれる. しかし, 温度だけでは回転子の高速駆動は不可能である. そこで我々は, 固定子の設計による遷移状態の安定化に加えて, 電磁波などの外部摂動によっ回転子て回転子を直接駆動すれば, 有効反応座標に沿った内部回固定子転運動を誘起して遷移状態を高速に通過させることができるとの着想に至った. この方針の下, 瀬高らはフェニレン
X = H or F
si = SiMe2環の片側をフッ素で置換した分子ジャイロスコープ（図 2
図 2 ジシロキサンかご状骨格をで X = F としたもの, 以下 RotF と表記）の合成に成功した持つ結晶性分子ジャイロスコー[7]. フッ素置換によって回転子−固定子間の立体反発が大プ. 無置換体（X = H）[4]を RotH,きくなり, RotH と比べて回転の活性化エネルギーは増加す
フッ素置換体（X = F）[7]を RotFると予測される. しかし, それと引き替えに RotF の回転子と呼ぶ.には大きな電気双極子モーメントが発生するため, 外部電場（特にテラヘルツ波）による回転駆動が原理的に可能である.
テラヘルツ波を利用した高速内部回転駆動の将来的な実現に向けて, まずは RotF の活性化エネ
ルギーの値や回転の動力学的機構を理論計算で解明する必要がある. ナノサイズの分子が集合した結晶性分子ジャイロスコープのシミュレーションを実行するためには, 以下に挙げる 3 つの条件を満たす計算法が必要である.
(ⅰ) 単位格子を無限に繰り返す周期境界条件の下での計算が行える.
(ⅱ) 計算負荷が低く, ピコ秒からナノ秒に及ぶ分子動力学シミュレーションが実行できる.
(ⅲ) 回転子−固定子間および分子間に働くファンデルワールス力までも高い精度で評価できる.
我々はこれら全ての条件を満たす密度汎関数強束縛（Density-functional-based tight-binding; DFTB）
法[8-11]を採用した. 固体物理学の分野で開発された DFTB 法は, 密度汎関数（Density functionaltheory; DFT）法に基づく半経験的手法である. 上述のように, 我々は既に DFTB 法を RotH のシミュレーションに適用してその有用性を確認している[5]. 本研究では, DFTB 法の代表的な計算プ
ログラム DFTB+[12]を東北大学サイバーサイエンスセンター所有の並列コンピュータシステム
NEC LX 406Re-2 に実装した. その並列化性能をテストした後に, RotF の結晶構造と回転動力学のシミュレーションを行った.

2. 密度汎関数強束縛（DFTB）法
最初に, DFTB 法の理論的枠組みを概説する. DFTB 法では, 系の電子密度 ρ を中性原子で構成される電荷の偏りが無い参照系の電子密度 ρ0 と電子密度揺らぎ δρ の和で表す. そして, ρ の汎関数である DFT 法の Kohn-Sham エネルギーを ρ0 の近傍で δρ の多項式として展開する. 1 次の項は 0となり, 系のエネルギーは参照系の Kohn-Sham ハミルトニアンと原子間反発から得られる 0 次の主要項, および δρ に依存する 2 次以上の補正項に分けられる.
電荷の偏りを無視して 0 次解を評価する手法を Non-self-consistent-charge (Non-SCC) DFTB 法[8,9]と呼ぶ. 強束縛近似（化学の用語では LCAO 近似）に基づいて参照系の Kohn-Sham 軌道を原子軌道の線形結合で表すと, Kohn-Sham 方程式は次の強束縛（永年）方程式に帰着される.

det  H 0   S   0(1)

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション— 3 —0H 0   H  は原子軌道間の Kohn-Sham ハミルトニアン行列, S   S  は重なり積分（内積）行列である. DFTB 法は, これらの行列要素を事前に DFT 法で見積もられたパラメーターとすることで, 高速かつ DFT 法に近い精度で結晶の計算を可能にする. (1)式を解いて Kohn-Sham 軌道の軌道
エネルギー  i  と展開係数 ci  が得られ, 系のエネルギーE はE  2 f i  i  Erep(2)iと与えられる. 2fi は Kohn-Sham 軌道 i の平均占有数であり, 通常 fi   0,1 は Fermi-Dirac 分布関数から決定される. Erep は原子間反発エネルギーであり, 解析的 2 体ポテンシャルの和として記述され, DFT 計算に合うようにフィッティングされる. 電荷の偏りを無視するこの Non-SCC DFTB 法
はナノカーボンなどの等核分子や非極性分子の計算に適している.
一方, 系のエネルギーに対して 2 次の補正項も考慮する手法を Self-consistent-charge (SCC)
DFTB 法[10]と呼ぶ. SCC DFTB 法では, (1)式の代わりに次の行列式を解く.

det  H   S   0(3)

H   H   は参照系の Kohn-Sham ハミルトニアン行列 H0 に各原子の電子密度揺らぎの間に働くCoulomb 相互作用を加えたものである.

10H   H  S    IK   JK  qK2K(4)

添字 I と J はそれぞれ原子軌道 μ と ν が局在する原子を指す. ΔqK は原子 K の Mulliken 電荷[13]の揺らぎ（中性原子との差）, γIK や γJK は原子の DFT 計算から決まる Hubbard パラメーター[14]と原子間距離に依存する関数である. (3)式と(4)式を自己無撞着的に解くためには, 反復計算が必要である. まず{ΔqK}の初期値を設定して(4)式の H を評価し, それを(3)式に代入して解く. 得られた展開係数 ci から新しい{ΔqK}を算出し, その値が収束するまで(3)式と(4)式を解く過程を繰り返す. 最終的な軌道エネルギー  i  を(2)式に代入して系のエネルギーE が求められる. SCC DFTB 法
は主に反復計算の数だけ Non-SCC DFTB 法よりも計算負荷が高く, 典型的には 10 倍以上に達する.
それでもパラメーターを使用しない DFT 法と比べればはるかに低負荷であり, 異核分子や極性分子の電荷分布を適切に記述することが期待できる.

 

3. MPI 並列版 DFTB+の並列化性能
DFTB 法に基づく量子化学計算パッケージの 1 つである DFTB+[12]は, ナノ材料などの大規模系への応用に向けた並列化が進められている. 本章では, 最近テスト版が公開された MPI 並列版
DFTB+の並列化性能を検証するために SiC 単結晶と RotF を対象として行った計算の結果を報告する. その準備として, 次節では RotF の結晶構造について説明する.

3.1

RotF の結晶構造

瀬高らが X 線結晶構造解析を用いて温度 273 K で観測した RotF の単位格子を図 3(a)に描画した. これは 4 分子で構成される直方晶（斜方晶）であり, 互いに直交する 3 つの格子ベクトルの長さは図中に記した通りである. 空間群 Pbcn に属し, 格子内のいずれか 1 つの分子に反転, 21 らせん, 対角映進の 3 種類の対称操作を施すことで他の 3 つの分子が得られる. 図 3(b)は格子内の 1分子を抜き出したもの, 図 3(c)はそれを固定子と回転子をつなぐ Si−C 結合（回転軸）に沿って眺めたものである. 一見すると, フェニレン環に 4 個のフッ素が置換され, 分子全体として C2 対称性を持っているように見える. RotF には図 3(d)に示した互いに完全に等価な安定構造 A と B が存

— 4 —

SENAC Vol. 48, No. 4（2015. 10）

在し, 結晶中の各分子が 1:1 の存在比でランダムに A または B の構造にある. X 線結晶構造解析では, 平均的な構造として 2 つの安定構造が重なって観測されている. 1 分子に含まれる原子数は195 個, 単位格子の 4 分子では 780 個に達する.
固定子を成す 3 本の鎖に特徴的な Si−O−Si（ジシロキサン）結合はほぼ直線的である. 図 3(c)の左右の鎖は捻れて広がっているのに対し, 中央の鎖は C2 軸方向へ伸びている. 一方, 回転子に着目すると, 環とフッ素の C−F 結合は環平面に対してわずかに C2 軸方向へ傾く.

図 3 X 線結晶構造解析から得られた温度 273 K における RotF の結晶構造[7]. (a) 4 分子で構成される単位格子. (b) 格子内の 1 分子. (c) 回転軸に沿って眺めた 1 分子. (d) 等価な 2 つの安定構造 Aと B の回転子近傍の拡大図. 観測された結晶構造は安定構造 A と B が重なった平均的な構造である.

3.2

SiC 単結晶と RotF を対象としたテスト計算

並列コンピュータシステム NEC LX 406Re-2 に MPI 並列版 DFTB+を実装した. コンパイルにはIntel Fortran Compiler XE 14.0 を使った. また, MPI ライブラリとして Intel MPI Library 4.1 Update 3
Build 048 を用い, MPI 並列用の線形代数ライブラリには Intel MKL Version 11.1 Update 1 に収録されている ScaLAPACK を指定した.
まず, SiC 単結晶をテストの対象とした. 原子数 1,000（原子価軌道数 No = 4,000）, 2,000（No =8,000）, 4,000（No = 16,000）の 3 つの系に周期境界条件を課して, SCC DFTB 法で電子エネルギー計算を行った. 計算終了までの経過時間（Wall clock time）とコア数の両対数グラフを図 4 に示した. コア数の増加と共に経過時間が期待通りにほぼ線形減少することが分かる.
次に, RotF に関しても経過時間のコア数依存性を調べた. 図 3 の結晶構造を基にして, 4 分子が— 5 —結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション

同じ安定構造にある単位格子を用意した. 図 5(a)は 4分子が全て安定構造 A にある場合, 図 5(b)は全て安定構造 B にある場合の単位格子を明示したものである.
両者は等価な結晶を成すので, どちらを計算に用いても同じ結果が得られる. 2×2×2 = 8 個の単位格子の集まりに対して周期境界条件を適用した. 総原子数は6,240, 総軌道数は 13,632 である. SiC 単結晶の 16,000軌道の場合と同様に, コア数の増加につれて経過時間が減少した（図 4 の○印）.
以上のように, 16,000 軌道程度の計算では, 384 コア（16 ノード）まで線形スケールすることが確認された.
しかし, 384 コアを使用した場合の並列化効率はわず
か 19 %であった（48 コアの場合は 50 %）. この原因を究明するために, 大規模計算パフォーマンス解析ツール Score-P[15]および Scalasca[16]を用いてボトルネック解析を行った. その結果, ScaLAPACK 収録の分割統治法（Divide-and-conquer algorithm）に基づくハミルトニアン行列の対角化サブルーチンにおいてデータの通信・同期を行う際に遅延が発生し, 計算負荷が使用されているコアに対して均等に分散されていないことが判明した. この通信･同期の遅延は, 子プロセス間のデータ送受信の順番が効率的に行われていないことが原因と考えられる. 現在, 子プロセス間のデータ送受信等を中心としたプログラムコードの最適化を進めている.

図 4 NEC LX 406Re-2 において測定した, 周期境界条件を課した SiC 単結晶の電子エネルギー計算の経過時間（Wallclock time）とコア数の関係 . 原子数1,000（原子価軌道数 No = 4,000）, 2,000（No = 8,000）, 4,000（No = 16,000）の結果をそれぞれ△, □, ◇で示した. ○はRotF の 2×2×2 = 8 個の単位格子の集まりに周期境界条件を適用した例. 実線はプロセス数を n 倍にしたときに経過時間が 1/n になる理想値を表している.

図 5 4 分子が同じ安定構造にある単位格子. (a) 全て安定構造 A の場合. (b) 全て安定構造 B の場
合. 回転子の向きが分かり易いように, 各分子のフッ素を点線で囲んだ. (a)を単位格子とする結晶と(b)を単位格子とする結晶は等価である.

4. RotF の結晶構造と回転動力学のシミュレーション
フッ素置換の影響で活性化エネルギーがどの程度増加するかを評価するために, RotF の有効反応座標に沿ったポテンシャルエネルギー曲線を求めた. また, 熱運動で環が回転して安定構造間— 6 —

SENAC Vol. 48, No. 4（2015. 10）を移動する動的挙動を探る分子動力学シミュレーションを実行した. 以下では, 比較的少数のコアで済む（並列化効率が高い）単位格子 1 個に周期境界条件を課した計算の結果を報告する.

4.1

結晶構造最適化

図 5 の 4 分子が同じ安定構造（A または B）にある単位格子 1 個に周期境界条件を課し, 格子定数は実験値に固定して SCC DFTB 法で構造最適化を行った. 図 3 の安定構造 A と B に対応する等価な 2 つの SCC DFTB 最適構造を図 6 に示した. 図 3 と比較すると, 固定子は C2 対称性を失い, Si−O−Si 結合が折れ曲がっている. 回転子では C2 軸に垂直だった環平面が傾いて C−F 結合とほぼ平行になった.

図 6 4 分子が同じ安定構造にある単位格子を SCC DFTB 法で最適化した結果. 格子内の 1 分子を抜き出した. (a)は安定構造 A, (b)は安定構造 B に対応しており, 両者は等価である. 紫色の 4 原子で回転軸周りの二面角を定義した.
回転子の配向を定量化する尺度として, 図 6 の紫色で指定した 4 原子で回転軸周りの二面角を定義し, その値を表 1 にまとめた. 構造 A と B は等価であるが, 前者の二面角は実験値に近いのに対し, 後者はややその差が大きい. 鎖の変形と環平面の傾きの影響が, A では相殺する一方, B では共に二面角を小さくするように働いているためである. 細部に関しては精度改善の余地があるが,分子全体の概形は図 3 と良く一致している.

4.2

表 1 安定構造における回転軸周りの二面角（deg.）.
AB

実験−76.9101.0

SCC DFTB−66.373.0

内部回転の有効ポテンシャルエネルギー曲線

上述のように, RotF の単位格子は 780 個の原子を内包する超多次元系である. 回転の活性化エネルギーを見積もるためには, その中から有効な反応座標を抽出する必要がある. そこで, 4 分子が同じ安定構造にある単位格子を起点として, 任意の 1 分子の二面角を変数として 1°ずつ変化させ, それぞれの角度において他の自由度を最適化した. 360 回の構造最適化には時間を要するので,この計算には Non-SCC DFTB 法を用いた（極小点と遷移状態のみ SCC DFTB 法でも求めた）.

— 7 —

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション

Potential energy (kcal/ml)

図 7 は二面角の関数として得られた内部回転の有効ポBA
テンシャルエネルギー曲線である. 80°の極小点では 4 分子が全て安定構造 B にあり, −65°の極小点は 1 分子のみ安定構造 A に移ったことに対応する. 環が伸びた鎖に重なる 0°付近と±180°付近に遷移状態が存在し, 活性化エネル
ギーの値はそれぞれ 5.6 kcal/mol と 4.5 kcal/mol（SCC DFTB法では 7.0 kcal/mol と 5.8 kcal/mol）である。フッ素で置換されたために RotH の値（0.7 kcal/mol）と比べて増加したが, それでも図 1 のかご状骨格を持たない分子ジャイロス
コープの中で最も活性化エネルギーが低いものと同程度
Dihedral angle (deg.)である.
図 7 内部回転の有効ポテンシャ回転に伴う固定子の構造変化を詳しく調べたところ,ルエネルギー曲線. 回転を起こす遷移状態の周辺では環との立体反発を避けるように鎖が分子が 2 つの極小点で成す構造を変形する. 特に, フッ素と鎖が接近する 0°付近では大きく明記した. 他の分子の構造は B.
歪む. RotF の有効反応座標は環が結合軸周りに回転するだけの単純なものではなく, 回転子と固定子が互いに相関した複雑な運動である. 一方, 2 つの極小点のエネルギーがほぼ等しい（SCC DFTB 法でも差は 0.4 kcal/mol）ことから, 分子間の双極子−双極子相互作用は小さいと示唆される. これは X 線結晶構造解析において回転子の配向秩序が現れない（ランダムに A または B の構造にある）事実と合致する. しかし, 反応座標の解析からは, 回転に際して固定子の変形が隣接する分子の構造にわずかながら影響することが分かった. 図 7 のポテンシャルエネルギー曲線が回転方向に対して非対称であることも結晶中の分子間相関を反映している.

4.3

熱運動の分子動力学シミュレーション

有限温度における熱運動が誘起する内部回転の挙動を解明するために, 全エネルギーが一定（小正準集団）の条件下で Non-SCC DFTB 法を用いた古典トラジェクトリ計算を実行した. 始時
刻の単位格子は 4 分子が全て安定構造 B にあると設定した. 運動温度の時間平均として系の温度を定義し, 300 K から 600 K の範囲を調査した.
300 K の場合, ピコ秒の時間領域では回転に至らなかった. 一方, 400 K では環がどちらの遷移状態も越えて回転することが確認された. 図 8(a)は 400 K の典型的なトラジェクトリにおける回転軸周りの二面角の時間変化をプロットしたものである. 環が 200 ps の間に 2 つの安定構造間を移動しながら回転することが読み取れる. また, 古典トラジェクトリ計算においても, 環が遷移状態を越える瞬間に固定子が変形して回転を促す様子が見られた.

図 8 古典トラジェクトリ計算の結果. (a) 温度 400 K の典型的なトラジェクトリにおける回転軸周りの二面角の時間変化. 有効ポテンシャルエネルギー曲線の極小点の位置を縦線で示した. (b)±180°付近の低い遷移状態を越える移動に関するアレニウスプロット.

— 8 —

SENAC Vol. 48, No. 4（2015. 10）環が安定構造に滞在する時間と回転に要する所要時間を比較すると, 前者の方が十分に長いとみなせる. そこで, 安定構造間の移動の速度定数 k を平均滞在時間の逆数として定義し, アレニウスプロットを作成した. 図 8(b)は±180°付近の低い遷移状態を越える移動に関するアレニウスプロットである. 直線でフィッティングして得られた活性化エネルギーは約 4 kcal/mol であり, 図 7 の有効ポテンシャルエネルギー曲線から算出した値に近い. この直線を外挿することにより, 室温（298 K）における速度定数を k = 4×108 s−1 程度と予測した. これは約 3 ns に 1 回の頻度で低い遷移状態を越えて 2 つの安定構造間を移動することを意味する. 実験では, RotF の内部回転は室温で NMR の時間スケールよりはるかに速いことが分かっており, 回転子−固定子間の立体反発が著しく大きい塩素置換体 RotCl では 83 μs に 1 回の頻度と見積もられている[7]. 我々の計算はこれらの実験結果と矛盾しないものである.

5. まとめと展望
次世代の光機能性材料として有望視される結晶性分子ジャイロスコープ RotF の回転動力学を解明するために, DFTB 法によるシミュレーションを行った. 計算結果の詳細な解析から得られた成果を以下にまとめる.
(ⅰ) 周期境界条件を課した DFTB 法は RotF の結晶構造を半定量的に再現できた.
(ⅱ) 内部回転の有効ポテンシャルエネルギー曲線の探索から, 2 つの遷移状態（活性化エネルギーは 5.6 kcal/mol と 4.5 kcal/mol）の存在が明らかになった.
(ⅲ) 環の回転に際しては固定子も相関して動く.
(ⅳ) アレニウスプロットから室温での回転頻度を算出し, 約 3 ns に 1 回の頻度で低い遷移状態を超えると予測した.
今後は, テラヘルツ波照射による RotF の回転駆動シミュレーションを行う計画である. RotH と比較することで, テラヘルツ波から回転子へのエネルギー供給に対するフッ素置換の効果を調べていく. 分子ジャイロスコープは分子機械としても位置づけることができ[17]，今後のさらなる展開が期待されている. このような大型計算機を活用した実在系のシミュレーション研究から, 高い機能性を有する分子の設計に繋げていきたい.

謝辞
本研究の計算の一部は, 東北大学サイバーサイエンスセンターの並列コンピュータシステム
NEC LX 406Re-2 を利用することで実現できた. なお, 本研究の一部は JSPS 科研費 15KT0138 および文部科学省 HPCI 戦略プログラム分野 2「新物質・エネルギー創成」の助成を受けて行われた.

参考文献[1] Z. Dominguez, H. Dang, M. J. Strouse, and M. A. Garcia-Garibay, Molecular “Compasses” and“Gyroscopes.” III. Dynamics of a Phenylene Rotor and Clathrated Benzene in a Slipping-Gear Crystal
Lattice, J. Am. Chem. Soc. 124, 7719 (2002).
[2] T.-A. V. Khuong, G. Zepeda, R. Ruiz, S. I. Khan, and M. A. Garcia-Garibay, Molecular Compasses
and Gyroscopes: Engineering Molecular Crystals with Fast Internal Rotation, Cryst. Growth Des. 4,15 (2004).
[3] Z. J. O’Brien, A. Natarajan, S. I. Khan, and M. A. Garcia-Garibay, Synthesis and Solid-State
Rotational Dynamics of Molecular Gyroscopes with a Robust and Low Density Structure Built with aPhenylene Rotator and a Tri(meta-terphenyl)methyl Stator, Cryst. Growth Des. 11, 2654 (2010).
[4] W. Setaka, S. Ohmizu, C. Kabuto, and M. Kira, A Molecular Gyroscope Having Phenylene Rotator
Encased in Three-spoke Silicon-based Stator, Chem. Lett. 36, 1076 (2007).

結晶性分子ジャイロスコープの構造と回転動力学の密度汎関数強束縛法によるシミュレーション— 9 —[5] A. B. Marahatta, M. Kanno, K. Hoki, W. Setaka, S. Irle, and H. Kono, Theretical Investigation of theStructures and Dynamics of Crystalline Molecular Gyroscopes, J. Phys. Chem. C 116, 24845 (2012).
[6] W. Setaka and K. Yamaguchi, Thermal modulation of birefringence observed in a crystalline
molecular gyrotop, Proc. Natl. Acad. Sci. USA, 109, 9271 (2012).
[7] W. Setaka, S. Ohmizu, and M. Kira, Molecular Gyroscope Having a Halogen-substituted p-Phenylene
Rotator and Silaalkane Chain Stators, Chem. Lett. 39, 468 (2010).
[8] D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner, Construction oftight-binding-like potentials on the basis of density-functional theory: Application to carbon, Phys.
Rev. B 51, 12947 (1995).
[9] G. Seifert, D. Porezag, and Th. Frauenheim, Calculations of Molecules, Clusters, and Solids with a
Simplified LCAO-DFT-LDA Scheme, Int. J. Quantum Chem. 58, 185 (1996).
[10] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, and G. Seifert,
Self-consistent-charge density-functional tight-binding method for simulations of complex materials
properties, Phys. Rev. B 58, 7260 (1998).
[11] M. Elstner and G. Seifert, Density functional tight binding, Phil. Trans. R. Soc. A 372, 20120483(2014).
[12] B. Aradi, B. Hourahine, and Th. Frauenheim, DFTB+, a Sparse Matrix-Based Implementation of theDFTB Method, J. Phys. Chem. A 111, 5678 (2007).
[13] R. S. Mulliken, Electronic Population Analysis on LCAO–MO Molecular Wave Functions. I, J. Chem.
Phys. 23, 1833 (1955).
[14] J. Hubbard, Electron Correlations in Narrow Energy Bands, Proc. R. Soc. A 276, 238 (1963).
[15] A. Knüpfer, C. Rössel, D. an Mey, S. Biersdorff, K. Diethelm, D. Eschweiler, M. Geimer, M. Gerndt,
D. Lorenz, A. Malony, W. E. Nagel, Y. Oleynik, P. Philippen, P. Saviankou, D. Schmidl, S. Shende, R.
Tschüter, M. Wagner, B. Wesarg, and F. Wolf, Score-P: A Joint Performance Measurement Run-Time
Infrastructure for Periscope, Scalasca, TAU, and Vampir, in Tools for High Performance Computing
2011, 79 (Springer, 2012).
[16] M. Geimer, F. Wolf, B. J. N. Wylie, E. Ábrahám, D. Becker, and B. Mohr, The Scalasca performancetoolset architecture, Concurrency Computat.: Pract. Exper. 22, 702 (2010).
[17] S. Erbas-Cakmak, D. A. Leigh, C. T. McTernan, and A. L. Nussbaumer, Artificial Molecular Machines,
Chem. Rev. 115, 10081 (2015).