SENAC Vol. 45, No. 3（2012.7）— 9 —［共同研究成果］

����������������������������� MPI �������
岩長 祐伸
物質・材料研究機構，科学技術振興機構さきがけ
周期長が光の波長程度である人工周期構造体（フォトニック結晶）に構成要素として金属ナノ構造を含む一般的なフォトニック結晶における光学応答を高精度に数値計算するためには，巨大なメモリを使用しながら 10000×10000 程度の一般複素数値行列の演算を大量に実行する必要がある．数年来，この計算コードの高速化，並列化に漸進的な改良を行ってきたが，さらなる高速化を実現するために MPI 化を実施した．本稿では，MPI 化によって得られた高速化の現状を中心に報告する．

1.

����

電磁波の数値計算に関しては多くのソフトウェアが存在し，商用シェアウェアで利用できるものとしてマクスウェル方程式を時間差分して電磁場を計算する FDTD（Finite Difference TimeDomain）法，空間を有限要素で分割してマクスウェル方程式を解く有限要素法，マクスウェル方程式をフーリエ変換した方程式を解いて反射率・透過率などを算出する RCWA（RigorouslyCoupled Wave Approximation）法などが比較的よく知られている．
いずれの方法でもマクスウェル方程式を数値的に解くことになるので，空間のグリッド分割が必要になる．定性的には，フォトニック結晶内における電磁波の波長 λ の数十分の１程度のグリッドの細かさが必要となる（RCWA 法では空間座量から電磁波の波数へのフーリエ変換が行われるので，さらにフーリエ波数 2πim/λ も大きな m 次まで取ることになる）．フォトニック結晶の単位胞に金属ナノ構造が含まれていると，ナノ構造のサイズ程度で電磁波が固有の分布を形成するため，ナノ構造の数十分の１程度のグリッドというさらに細かい空間分割が必要となってくる．
図１は金属フォトニック結晶の細線で外枠を示す単位ドメイン内における金属部表面を分割した一例を示している．この例は有限要素法によるもので境界 yz，zx 面に周期境界条件を課して計算を行う．円柱状のロッドは直径 150 nm，x 軸方向の周期は 900 nm，z 方向の金属の厚さは 150 nmである．金属ナノロッドの直径 150 nm に対して，その数十分の１は 2~5 nm になる．

図１

一般フォトニック結晶単位ドメイン（細線）内の金属部分のグリッド分割の例

— 10 —

SENAC Vol. 45, No. 3（2012.7）

FDTD 法や RCWA 法においても単位ドメインに対応する分割を設定して，対象の光学応答を算出することになる．結果として周期構造の単位ドメインを 108 程度の要素に分割して，各グリッドに物質パラメータである複素誘電率を割り当て，マクスウェル方程式の膨大な数値計算（最終的には数値行列演算）を行うことになる．
研究の現場で未知のフォトニック結晶の性質を明らかにするためには，試行錯誤を含めて多くの数値計算を実行しながら，計算精度を追求することが必要となるので，筆者は RCWA 法の原著論文［1］を基に有限次切断フーリエ展開を高速に取り扱うアルゴリズムを取り入れ，周期構造が積層した 3 次元的なフォトニック結晶を数値的に安定に取り扱うことができる散乱行列法［2］を組み込んだ数値計算コードを開発し，SX-9 上で運用してきた．これまでの RCWA・散乱行列コードのアルゴリズムの詳細，運用実例，改良については本誌上で報告してきた［3－7］ので繰り返さないが，漸進的に高速化・並列化の向上を実現し，MPI 化を除けば，改良の上限に達していた．
計算機環境についても読者の関心があると思われるので，ここで簡単に触れておきたい．近年のワークステーションの進歩に伴い，100 GB 程度のメモリをもつワークステーションを研究室に備えて計算を実行していくことも選択肢としてありうるが，RCWA 計算では対象次第では 1 TBに近いメモリを要する場合も稀ではないため SX-9 上でのコード運用に大きな利点を筆者は感じている．また，メモリだけならばクラスター型の計算機なら対応できるが，数値行列演算に関して SX-9 のほうが数倍から数十倍速かったという経験もある．コード自体に多くのチューニングを施さないと高速化できないクラスター計算機はユーザーにとっては恩恵を感じにくいと想像する．一方，SX-9 上では自動並列化オプション（-Pauto）のみで相当程度の並列化率を得て高速化できる利点がある［7］．
本稿では，以上の状況を受けて，今回 RCWA・散乱行列コードの MPI 化とその結果について述べる．通常，フォトニック結晶というと，シリコンなどの半導体材料または透明誘電体であるオパールなどからなるものを通常想定するので，以下では金属ナノ構造も含む，より一般的なフォトニック結晶を単に一般フォトニック結晶と呼ぶことにする．

2．MPI ����一般フォトニック結晶の光学応答では線形過程によって生じるもののみを考える．この場合，反射光，透過光に加えて回折光が生じうる．これらの応答の特性はスペクトルとして表現し，その性質を吟味することになる．図２は図１で示した金属フォトニック結晶の模式図と反射（R）スペクトル（赤線），透過（T）スペクトル（青線）である．入射角度 θ が 0 度（実線）と 10 度（点線）の場合を示している．

図２ 金属フォトニック結晶の模式図と反射（赤線）・透過（青線）スペクトル— 11 —

金属ナノ構造を含む一般フォトニック結晶の光学応答計算コード MPI 化による高速化

図２左では入射光の配置を示している．入射面（入射光と反射光の進行ベクトルが張る平面）
は xz 面であるように設定し，入射偏光は電場ベクトル Ein が xz 面内にある p 偏光とした．対応する実験も行い，計算結果とよい一致を示すことを明らかにした［8］．図２右では横軸を光の波長で表示した反射・透過スペクトルを示し，入射角度によってピークやディップが大きく変化していることが分かる．これらのスペクトル特性から共鳴状態を明らかにするのが筆者の数値計算を通じた物理学的な主題であり，実際に起源の異なるプラズモン共鳴として分類できるのであるが，考察結果と構造パラメータの詳細は文献［8］に譲り，ここでは数値計算の条件などについて述べていく．
図２の反射・透過スペクトルは RCWA 法によって SX-9 上で計算した結果であり，xy 面を 5×5nm2 のグリッドで分割して，フーリエ波数を±20 次まで取った．この計算を p16（16 CPU 下）で実行した結果，使用したメモリは 59.6 GB でスペクトル上の１点を計算するのに 9.8 分要した．
１つのスペクトルには 139 点のデータ点があるから，22.7 時間かかったことになる．
光学スペクトルの各点は線形過程では互いに独立であるから，１つのスペクトル（列データ）を計算するために MPI 化で強制的に並列化することは大幅な計算時間の短縮につながることが期待できる．単純に言えば，図２の例を p64 で 16 CPU×4 の実行を行うことができれば，計算時間は 1/4 に短縮できる．この大幅な時間短縮，つまり高速化を実現するために今回コードの MPI化を実施した．
MPI 化の方針は図３に図式化した．基本的な方針としては，スペクトル計算の主要ループを強制的に並列化すればよい．ただし，主要ループの前で共通保持できるポインターは確保し，物質パラメータなどを割り当てている．

図３

MPI 化の概念図

計算実行時のメインループは図３で簡略化して示しているように E-LOOP が担う．このループはスペクトルの横軸に相当する量を掃引する．主な計算処理は，
1. フーリエ変換したマクスウェル方程式の係数行列を構成する CONSTRUCT_MATRIX

— 12 —

SENAC Vol. 45, No. 3（2012.7）2.
3.

各周期層のなかのマクスウェル方程式を解く SOLVE_MATRIX
各周期層の固有モードを使って一般フォトニック結晶全体の応答を表現する散乱行列を構成する CONSTRUCT_S_MATRIX
4. 最後にインプット（入射光）に対してアウトプット（反射光，透過光など）を散乱行列から算出する CAL_R_Tがある．各処理で使う配列が E-LOOP ごとに独立であるから，今回の MPI 化は実行できた．物理量として対象の線形応答からスペクトルを算出する場合は，今回のように電磁波の問題に限らず，同様の MPI 化が一般に可能である．
図３の左側は自動並列化のみでの実行を想定しており，p16 または p8 または s での実行に対応する．一方で図３右は MPI 化後に p64 に 16 CPU×4 並列で実行することを模式化している．スペクトルのデータ点が 8 点であれば，2 点×4 並列となり，理想的には計算時間は 1/4 になる．一般にはデータ点が 4 で割り切れるとは限らないが，その場合は端数を E-LOOP1 と E-LOOP2 に割り振る実装になっている．スペクトル計算において通例データ点は 100 から 1000 点ほどあるので，実際の運用において端数が出たことによる遅延が問題になることはなかった．
表１は p16 と p32（16 CPU×2 並列），p64（16 CPU×4 並列）での計算時間と使用メモリの変化の典型的な例である．テストケースとして 3 層積層からなる一般フォトニック結晶の反射・透過スペクトルを計算した（図１，２の構造とは異なる）．
p16
p32
13.2
3.38
データ１点の計算時間（分）
50.15
100.67
使用メモリ（GB）
表１ MPI 化の有無による計算時間，使用メモリの変化

p64
1.78
200.34

表１から，p64 での計算時間が 1.78/13.2 = 1/7.4 と想定の 1/4 よりさらに高速化しているように読み取れる．この例では p16 でのデータ点数を 21 点に抑えたため，偶然時間のかかる波長域で計算してしまった可能性がある．共鳴状態の波長においては SOLVE_MATRIX のなかで call している行列の固有値・固有ベクトルを求める関数の実行に時間がかかる傾向がある．表１以外の例では，およそ 1/4 程度の時間短縮が得られており，当初の目標通りの高速化が実現できていると言える．また，p32 と p64 を比較すると，1.78/3.38 = 1/1.9 と p64 で約 1/2 の実行時間で計算を実行できている．
使用メモリに関しては，MPI による並列化で同じ配列を並列数だけ確保するので，p16 に対して p64 では約 4 倍，p32 では約 2 倍のメモリを使用している．これも想定通りであり，MPI の実装がうまくできていることの証左である．

3. メタ��リ��������の��
この節では SX-9 上での光学応答計算コード運用から得られた，この１年の結果［8－10］のなかからメタマテリアルに関する成果［9,10］について述べる．メタマテリアルとは，一般フォトニック結晶のなかで考察する波長域よりも周期長が小さいものであり，回折光が生じないことから，従来の固体媒体では実現困難な性質を備えた新しい電磁波伝播媒体となることが期待されている．

3.1 ���������メタ��リ��������モー�の分�計算
フィッシュネット・メタマテリアルとは，金属・絶縁体・金属の積層構造に貫通孔を周期的にあけたフォトニック結晶のことを指す．図４左は構造模式図であり，図４右は p 偏光下の波数・光エネルギー分散図を示している．分散のデータ点は光学応答から得られた吸収スペクトルのピーク位置をプロットすることで得た．光吸収量 A は入射光の電磁エネルギーを 1 と規格化したと— 13 —

金属ナノ構造を含む一般フォトニック結晶の光学応答計算コード MPI 化による高速化

きに A  1  R  T で表される（ R ：反射率， T ：透過率）．最低次下枝 a の分散式は

   k x  1(1)と表される．ただし，   0 ， k x は波数， 1 は k x  0 での最低次のエネルギーである．

図４ フィッシュネット・メタマテリアルにおける固有モードの波数・エネルギー分散
金属層に挟まれた絶縁体層内部の導波路モードの群速度 vg はvg k x(2)で与えられ，式（1）から v g    0 となる．つまり，最低次下枝 a の導波路モードは負の群速度を持つことが分かった．このような特異な条件が満たされることで，このフィッシュネット・メタマテリアルは，斜め入射下で負の屈折現象が生じる媒体となっていることが今回初めて明らかにできた［9］．なお，最低次のモードであっても上枝 b の群速度は正であり，負の屈折現象は生じない．

3.2� ������������������������������
金属ナノロッドが周期的に配列した場合に，単独のロッドだけでは生じない集団的な共鳴状態を見出した．図５はその具体的な構造と透過スペクトルを示している．周期構造に上から入射光が当たる配置を考える．
金属ナノロッド間の距離（ギャップ）を 0 nm から 20 nm まで少しずつ変えていく．ギャップが 0 nm のときは入射直線偏光が ψ=45 度（ψ は x 軸と直線偏光のなす角）に対して波長 1500 nm以上で透過率がほぼ 0 になっている．一方，入射偏光 ψ=135 度に対しては同じ波長域で 40%程度の透過率がある．このことは ψ=45 度では広帯域な共鳴状態が形成されて，光エネルギーが吸収されていることを示唆している．実際，電磁場分布の解析からナノロッド間の集団的な共鳴が生じていることを明らかにした［10］．ナノロッド間の集団共鳴を示したのは，筆者の知るかぎりにおいてこの例が初めてである．この集団共鳴はロッド間の距離に非常に敏感でわずか 10 nm のギ

— 14 —

SENAC Vol. 45, No. 3（2012.7）

ャップでもほぼ失われてしまう．構造設計の精度が非常に重要であることを示唆している．

z

図５ 金属ナノロッド列における構造変化と偏光透過スペクトル

金属ナノ構造を含む一般フォトニック結晶の光学応答計算コード MPI 化による高速化— 15 —

一方で入射偏光 ψ=135 度のとき，つまり入射直線偏光がナノロッド列と直交しているとき，透過スペクトルはほとんどギャップに依らない．この結果は集団的な共鳴状態が偏光選択性をもつことを示している．このようにフォトニック構造と光の状態（偏光）の組み合わせを含めて，対象の光機能を設計していく必要がある．

4. ࡲ࡜ࡵ࡜௒ᚋࡢᒎᮃ
一般フォトニック結晶の光学応答を数値計算するコードの MPI 化を実装し，高速化（計算時間の短縮）を試みた．当初の期待通り，p64 では p16 の約 1/4 の実行時間でコードの運用ができるようになった．光学スペクトルの計算が電磁波の線形応答をもとにしたものであったので，各データ点を独立に扱うことができ，自ずと MPI 化に適していた．電磁気学の問題に限らず，線形スペクトル計算の MPI 化は物理学一般において有効であると考えられる．
計算時間の短縮は様々な対象をより多く扱うことを可能にするので，コードを運用するユーザーとしても大きな恩恵にあずかることができる．この１年だけでも様々な対象を詳細に研究できた［8－10］ことも端的な成果と言える．
本稿で紹介した例では p16 実行下で 50 GB 程度のメモリを使ったものが多かったが，その他の事例ではメモリ 500 GB というものもあった．SX-9 でもメモリの上限をそろそろ意識しなくてはならない領域に踏み込みつつある．私見ながら，SX-9 のスペックは電磁気学の問題を解くのに適していると感じている．量子力学の問題の多くは現在の大型計算機を持ってしても解けていない
ことは周知であり，チャレンジングな課題であることは確かであるが，単に計算機のスペック不足なのか，それとも良いアルゴリズムを見つけることができていないだけなのか，筆者にはいまだに判然としない．今回用いた RCWA 法に関して，実用に耐えるアルゴリズムの発見［1］までに約 30 年を要した歴史的経緯がその思いを強めている．
今回の MPI 化の効果として，さらに膨大な計算量を必要とする数値的な新規フォトニック構造探索が可能になった．遺伝アルゴリズム［5］と今回の光学応答計算コードの融合による構造探索を現在実行している．比較的単純な構造に潜む共鳴状態を解くことがこれまでの主な研究であったが，人間の想像を超える新しい人工ナノ構造の発見につなげていくことを目標に SX-9 を今後活用していきたいと考えている．

ㅰ㎡
本研究における MPI 化の実装は SX-9 開発元 NEC のご協力を得て，東北大学サイバーサイエンスセンターと共同で実施したものであり，期待にたがわぬ MPI 化を行っていただきました．この場を借りて厚くお礼を申し上げます．また，本研究の一部は科学技術振興機構さきがけ，科学研究費補助金（No. 22760047）の支援を受けて行われました．
ཧ⪃ᩥ⊩
[1] L. Li, “New formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. Soc.
Am. A 14 (10), 2758–2767 (1997).
[2] L. Li, “Formation and comparison of two recursive matrix algorithm for modeling layered diffractiongratings,” J. Opt. Soc. Am. A 13 (5), 1024–1035(1996).
[3] 岩長祐伸,「散乱行列法を用いたフォトニック結晶の光学応答解析」SENAC 39 (3), 25–32(2006).
[4] 岩長祐伸,「メタマテリアルにおける有効光学定数の決定法と応用」SENAC 40 (3), 5–14 (2007).
[5] 岩長祐伸,「遺伝アルゴリズムを用いた光機能性人工構造体の探索」SENAC 41 (3), 43–51(2008).
[6] 岩長祐伸，「メゾ周期構造体における電磁波散乱の高精度数値計算」SENAC 42 (4), 9–18 (2009).
[7] 岩長祐伸，「積層プラズモニック結晶における光機能性発現」SENAC 44 (2), 49–56 (2011).
[8] M. Iwanaga, N. Ikeda, and Y. Sugimoto, “Enhancement of local electromagnetic fields in plasmonic

— 16 —SENAC Vol. 45, No. 3（2012.7）

crystals of coaxial metallic nanostructures,” Phys. Rev. B 85 (4), 045427 (2012).
[9] M. Iwanaga, “In-plane plasmonic modes of negative group velocity in perforated waveguides,” Opt.
Lett. 36 (13), 2504–2506 (2011).
[10] M. Iwanaga, “Collective Plasmonic States in Metallic Nanorod Array and Their Application,” inNanorods, edited by O. Yalcin (InTech, Rijeka, 2012) Chap. 4. (http://dx.doi.org/10.5772/36198)