東京工業大学 TSUBAME 共同利用
令和３年度利用終了課題 利用成果報告書集

東京工業大学 学術国際情報センター
共同利用推進室https://www.gsic.titech.ac.jp/tsubame本報告書集および個別の報告書の PDF ファイルは、以下の URL にあります。
令和３年度採択実績および利用終了課題報告書
https://www.gsic.titech.ac.jp/kyodou/kadai_r3■令和３年度産業利用 利用成果報告書 一覧
申請課題名
所属機関／利用課題責任者

頁

車載レーダにおけるターゲット表面電流分布の FDTD 法による数値解析
1
マツダ株式会社技術研究所／山本雅史
分子標的型農薬開発のためのタンパク質-阻害剤複合体シミュレーション
5
株式会社アグロデザイン・スタジオ／田中良樹

■令和３年度学術利用 利用成果報告書別紙の提出免除課題 一覧（利用成果を論文/学会等にて発表した要旨等の提出により提出免除）
申請課題名
所属機関／利用課題責任者
論文/学会等における発表済利用成果の情報
計算化学による固体分散体中の薬物とキャリアの相互作用の評価
千葉大学大学院薬学研究院製剤工学研究室／東顕二郎[1] Computational approach to elucidate the formation and stabilization mechanism of amorphous formulation
using molecular dynamics simulation and fragment molecular orbital calculation. Xiaohan Ma, Kenjirou
Higashi, Kaori Fukuzawa, Keisuke Ueda, Kazunori Kadota,Yuichi Tozuka, Etsuo Yonemochi, Kunikazu Moribe,International Journal of Pharmaceutics, (査読あり), 615, 121477 (2022)
DOI:10.1016/j.ijpharm.2022.121477
大規模メタサーフェスの光学応答解析
東京農工大学／岩見健太郎[1] Chikara Ogawa, Sotaro Nakamura, Takumi Aso, Satoshi Ikezawa, and Kentaro Iwami, “Rotational varifocal
moiré metalens made of single-crystal silicon meta-atoms for visible wavelengths,”Nanophonics, Accepted,
2022(査読あり) DOI: 10.1515/nanoph-2021-0690[2] Kentaro Iwami, “Rotational varifocal Moire metalens at the visible Wavelength”, The 5th A3 (China- JapanKorea) Metamaterials Symposium, Online, June 27th 2021
[3] Kentaro Iwami, “Dielectric metasurface for imaging and visualization”, International Symposium on Imaging,Sensing, and Optical Memory (ISOM) 2021, Online, October 6th 2021
MPI アプリケーションのプロファイルおよびトレース予測に関する研究
電気通信大学／三輪忍[1] 岡田 悠希，三輪 忍，八巻 隼人，本多 弘樹， MPI における小規模実行時の通信トレース解析による大規模実行時の通信タイミング予測の評価， 情報処理学会研究報告 2021-HPC-182, No.16, pp.1-8(2021)
第一原理計算による SOFC 電極上の炭素析出解析立命館大学理工学部／渡部弘達[1] First-Principles Calculation and Multi-Scale Observation for a Carbon Deposition on a SOFC Anode Near
the Interface, Hirotatsu Watanabe, Teppei Ogura and Katsunori Hanamura
DOI:10.1149/10301.0825ecst
分子動力学シミュレーションによる新興再興感染症の研究
国立感染症研究所病原体ゲノム解析研究センター／横山勝[1] Kotani O, Suzuki Y, Saito S, Ainai A, Ueno A, Hemmi T, Sano K, Tabata K, Yokoyama M, Suzuki T, Hasegawa
H,Sato H. Structure-Guided Creation of an Anti-HA Stalk Antibody F11 Derivative That Neutralizes BothF11-Sensitive and -Resistant Influenza A(H1N1)pdm09 Viruses. Viruses. 2021 Aug 31;13(9):1733.
DOI:10.3390/v13091733. PMID: 34578314; PMCID: PMC8473006.

申請課題名所属機関／利用課題責任者
論文/学会等における発表済利用成果の情報
多様なデータを活用する深層学習モデルの検証
豊橋技術科学大学／後藤仁志[1] 立花尚登、五十幡康弘、後藤仁志 （豊橋技科大），配座データを用いたニューラルネットワークによる分子物性予測，日本コンピュータ化学会 秋季年会 ポスター 3P08, 2022、長野■令和３年度学術利用 利用成果報告書 一覧
申請課題名
所属機関／利用課題責任者

頁

動的均衡にある神経回路の機械学習とその場の理論的解析
7
日本大学医学部薬理学分野／早川隆
大規模ソフトウェア分散共有メモリシステムにおける高性能計算の研究
11
成蹊大学理工学部情報科学科／緑川博子
廃棄物最終処分場における間隙内流体挙動の数値解析
17
埼玉県環境科学国際センター／鈴木和将
揺動電磁気シミュレーションを用いた光放射現象に関する研究21
新潟大学工学部工学科／櫻井篤
液体金属流れ CFD 手法の開発及び核融合研究への応用
23
九州大学応用力学研究所／胡長洪
海洋底探査を加速する AI とシミュレーション技術の開発
25
大阪府立大学大学院工学研究科／橋本博公
逐次導入法による長鎖化合物の全原子自由エネルギー計算
27
大阪大学大学院基礎工学研究科／松林伸幸
ヒト血清アルブミン-環状ペプチド複合体の相互作用解析
29
長岡工業高等専門学校／和久井直樹
電気コイル隙間内に流れ込む冷却液挙動に関する大規模数値解析
31
大阪府立大学工学研究科／金田昌之
ヘロナミド類の抗真菌作用メカニズムの解明のための計算分子設計技術の開発
33
北陸大学薬学部／齋藤大明
LRnLA アルゴリズムを用いた物理シミュレーション
35
法政大学情報科学部／善甫康成
GPU クラスタを用いたミリ波帯・テラヘルツ帯大規模広帯域電波伝搬シミュレーション
39
国立研究開発法人情報通信研究機構／チャカロタイ ジェドヴィスノプ
フェーズフィールド法に立脚した多結晶組織形成の大規模シミュレーション
43
東京農工大学大学院工学研究院／三好 英輔

申請課題名所属機関／利用課題責任者

頁

水圏動物が有する酵素の構造解析
47
東京大学大学院農学生命科学研究科／潮 秀樹
蛋白質天然変性領域の動力学特性に関する計算科学的検討
49
立命館大学生命科学部／笠原浩太
超流動ヘリウムにおける量子乱流の数値的研究
51
慶應義塾大学自然科学研究教育センター／湯井悟志
新規 BRCAness 誘導阻害剤開発に向けた理論研究
55
名古屋大学大学院医学系研究科／小関準
粗さを用いた熱・運動量輸送の非相似制御に関する研究
57
大阪府立大学大学院工学研究科／桑田祐丞
二酸化炭素の分離回収に適した分離材料の探索
61
産業技術総合研究所／藤井達也
光熱変換機構に関する第一原理計算
63
山形大学大学院理工学研究科／江目宏樹
シアシックニング機構の解明に向けた理論シミュレーション
65
高度情報科学技術研究機構／山田達矢
二酸化炭素の分離回収におけるゼオライト構造の最適化
71
産業技術総合研究所／池田歩
臨床情報統合データベースの機械学習解析73
日本大学医学部生体機能医学系薬理学分野／浅井聰
CFD 解析を用いた卓球ボールの空力特性の計算
75
新潟工科大学／伊藤建一

（様式 20）成果報告書

TSUBAME 共同利用 令和 3 年度 産業利用 成果報告書

車載レーダにおけるターゲット表面電流分布の FDTD 法による数値解析
Numerical analysis of surface current distribution for an automotive radar by FDTD method
山本 雅史
Masashi Yamamoto
マツダ株式会社
Mazda Motor Corporation
http://www.mazda.co.jp
レーダ電波が照射された車両の表面電流分布は、車両のどの箇所からレーダ電波が強く反射しているかを把握する上で重要な物理量である。表面電流分布の数値解析は、時間領域でマクスウェル方程式を解く FDTD 法（Finite
Difference Time Domain method）により可能である。しかし、車載レーダのミリ波・準ミリ波帯で車両を対象にした
FDTD 計算を実施するためには、大規模な計算リソースと数値分散誤差の問題があった。今回、TSUBAME3.0 とFDTD(2,4)法を活用することで問題解決でき、実物大のターゲット車両の表面電流分布を解析できた。これにより、車両からのレーダ反射波の発生メカニズムおよび受信電力の変動メカニズムを解明する手掛かりを得ることができた。
Surface current distribution is an important physical index for understanding from which part of atarget vehicle radar waves will be strongly reflected. Surface current can be calculated by FDTD method.
However, FDTD method requires large-scale calculation resources and has large numerical dispersion
errors in millimeter-wave and quasi-millimeter-wave bands analyses. These problems could be solved
by using TSUBAME3.0 and FDTD(2,4) method, to analyze the surface current distribution on a target
vehicle. This method will make us able to be obtained clues to elucidate mechanisms of reflected waves
generation and fluctuation of received power for a vehicle target.

Keywords: automotive radar, surface current, FDTD, electromagnetic simulations, vehicle target
1.

ヨー角変化に対して RCS 特性と受信電力が大きく変動

背景と目的
車載レーダは、自動運転システムや安全運転支援シすることが分かっている[2]。現状、ターゲット車両のヨ

ステムとって走行環境を認識するための重要なセンシ

ー角変化に伴う RCS 特性の変動メカニズムや車両特ングデバイスである。レーダ性能向上のためには、自有の反射波の発生メカニズム、受信電力の変動メカニ

車の先行車両を検知する場合や、自車の前方交差点

ズムは解明できていない。

を横断する車両を検知する場合など、あらゆる走行シ

電波がターゲットに入射すると、電磁誘導によりターーンを想定する必要がある。つまり、アンテナに対するゲット表面に電流が流れ、その電流により周囲に磁界ターゲット車両の相対角度（ヨー角）が 360 度変化すると電界が発生して反射波となる。したがって、ターゲットことを想定し、角度変化に伴う受信電力の変動を考慮

上の表面電流の生じ方を知ることができれば、照射さしてアンテナや信号処理などを開発する必要がある。

れたレーダ電波がどのようにターゲットで反射するのかRCS[1]はターゲット車両からのレーダ反射強度を把を把握できると考えられる。

握する上で重要な特性であるが、RCS の計算式ではタ

そこで我々は、ターゲット車両からのレーダ反射波のーゲットがアンテナに対して無限遠方に存在し、大きさ

発生メカニズムおよび受信電力の変動メカニズムを解の無い点と見なせることが前提である。車載レーダの明する手掛かりを得るために、車両の表面電流分布の場合、ターゲット車両までの距離は最大で 200m 程度

数値解析に取り組んだ。

であり、RCS 計算式の適用範囲外となる。

2.

車両の表面電流分布数値解析の課題

さ ら に 、 ミ リ 波 帯 ・ 準 ミ リ 波 帯 （ 周 波 数 ： 76GHz 、

車載レーダの電波伝搬シミュレーションは、幾何光学

24GHz）車載レーダの波長（4mm、12mm）に比べて、近似を用いたレイトレース法や PO 法、GO 法が一般的

ターゲット車両のサイズが大きく形状が複雑であるため、に用いられる。しかし、上記の数値解析手法は、電波

1

（様式 20）成果報告書の波動性を厳密に再現できないことと時間領域の解析た大規模クラスタ型スーパーコンピュータである。1 台

ではないことから、電磁誘導による表面電流を計算での計算ノードには、NVIDIA 製 GPU Tesla P100（単きない。FDTD 法（Finite Difference Time Domain

精度浮動小数点演算性能：9.3 TFLOPS、GPU メモ

method）[3]は、電波の支配方程式であるマクスウェル

リ：16GB）を 4 機搭載している。全計算ノードの GPU メ

方程式を差分化して、時間領域で電界と磁界を解く手

モリ量は約 34.6TB である。また、GPU P100（CUDA法である。この計算原理より、ターゲット表面の電界とコア数：3,584 コア、メモリ・バンド帯域幅：720 GB/s）は、磁界の時間変動から表面電流を計算できる。

ハイパフォーマンス CPU と RAM メモリ（CPU：Intel

しかし、FDTD 法でターゲット車両の表面電流を解

Xeon E5-2680 14 コア、RAM：DDR4-2400 メモリ帯

析する上で 2 つの問題点がある。1 つ目の問題点は、域幅 19.2GB/s）と比較して、数百倍の並列計算コア数

FDTD 法では、波長の 1/10 以下のサイズのセルで解と数十倍のメモリ帯域幅を有する。

析空間全体を分割して、時間ステップを進めながら繰り以上のような TSUBAME3.0 のマルチ計算ノードの返し計算するため、大規模な計算リソース（メモリ量、計

マルチ GPU を用いることで、FDTD 法の計算リソース

算時間）を必要とすることである。具体的には、ミリ波帯の問題を解決できた。今回取り組んだターゲット車両の（周波数：76GHz、セルサイズ：0.4mm）・準ミリ波帯

表面電流分布解析のための FDTD 計算は、約 20 億

（周波数：24GHz、セルサイズ：1.2mm）で車両を対象

セルの解析規模で 8 ノードの 32GPU を必要とした。計にした解析では、セル数が数十億セルから数千億セル算実施の結果、所要メモリ量は 458GB、計算時間は約となり、所要メモリ量は数百 GB から数 TB 程度、CPU

50 分であり、有効性を検証できた。

計算は数百時間を所要する見積りである。

4.

FDTD(2,4)法による位相誤差の解決

2 つ目の問題点は、FDTD 法が差分近似を用いた

一般的な FDTD 法は、マクスウェル方程式の時間微計算原理であるため、数値分散誤差（位相誤差）が生分と空間微分を 2 次精度で差分近似した Standard-

じることである。波長に対する伝搬距離が長くなるほど

FDTD 法である。この計算手法は、自由空間中で位相位相誤差は累積するため、高周波数帯の大規模空間

速度が光速とならず、伝搬距離が長くなるほど位相が解析では伝搬波形が崩れる。

遅れる。セルサイズをより小さくすることで位相誤差を

今回、上記 2 つの問題点に対して、TSUBAME3.0

低減できる。しかし、セルサイズを波長の 1/20 と設定しを活用し、我々が開発した高精度差分スキームのた場合で 1 波長分の距離の伝搬で約 1°の位相誤差

FDTD（2,4）法プログラムを用いてターゲット車両の表が生じ、50 波長分の伝搬では約 50°、100 波長分の面電流分布の数値解析に取り組んだ。

伝搬では約 100°の誤差が生じることとなり、許容でき

3.

TSUBAME3.0 による計算リソース問題の解決ない。またセルサイズをより小さくした場合、計算リソー

FDTD 法の計算アルゴリズムは、リープフロッグアル

スが指数関数的に増大するため計算リソースの問題をゴリズムで、時間ステップを進めながら 1 ステップ過去誘発する。

の電界と磁界から現在の電界と磁界を交互に計算する。

一方、FDTD(2,4)法は、時間微分を 2 次精度差分で、そして、空間領域を分割したセルに対して当該セルと隣空間微分を 4 次精度差分で近似する。この計算手法は、接したセルの 1 ステップ過去の電界と磁界を用いて、当時間領域で位相が進み、空間領域で位相が遅れる。こ該セルの電界と磁界を計算する。したがって、1 ステッの特徴から、時間領域の計算分解能（時間ステップ）に

プ過去の計算結果を保持することで空間領域全体を一作用するクーラン条件 CFL を最適化することで、時間括して並列計算ができる。ただし、セル毎の電界と磁界

領域と空間領域で生じる位相誤差を相殺できる。

データをメモリから読み込みと書き出しをするため、計

セルサイズを波長の 1/10 に設定した場合のクーラン

算速度はメモリ・バンド帯域幅（メモリ・データ転送速度）

条件 CFL と位相誤差の関係を図 1 に示す。この結果に律速される。

から、FDTD(2,4)法で CFL を 0.1363212 に設定する

TSUBAME3.0 は、540 台の計算ノードで構成されことで 1 波長分の伝搬の位相誤差が 1.09×10-5 °とな

2（様式 20）成果報告書
表 1. FDTD 法の計算条件
周波数
解析空間

10.0 [GHz]3.07 x 5.76 x 3.07 [m]（1024 x 1920 x 1024 [cell]）

セルサイズ

2.99x10-3[m] （≒1/10λ)

計算手法

FDTD(2,4)

図 1. 1 波長分の伝搬における位相誤差と CFL の関係
クーラン条件 CFL 0.1363212
時間ステップ1.36x10-12 [sec]放射源

平面波＋ガウスパルス波

計算ステップ数

7,200 [回]吸収境界PML 32 層、R0=1.0-32、M=4

変数の型

float (GPU)ターゲット表面電流は、時間領域でマクスウェル方程図 2. ターゲット車両の表面電流分布の解析条件式を解く FDTD 法で計算可能である。しかし、ミリ波帯・準ミリ波帯で車両を対象にした FDTD 解析を行うためり、ほぼ無視できるレベルに低減できた。これより、には、大規模な計算リソースが必要となることと、解析

FDTD 法の位相誤差の問題を解決できた。
5.

規模に応じて位相誤差が大きくなる問題があった。

ターゲット車両の表面電流分布の数値解析

実物大のターゲット車両を対象に、準ミリ波帯（周波上記の問題に対して、マルチ計算ノードのマルチ

数：10GHz、波長：30mm）レーダが照射されたときのGPU を搭載した TSUBAME3.0 を活用し、我々が開表面電流分布を FDTD(2,4)法で解析した。図 2 に、ヨ発した高精度差分スキームの FDTD（2,4）法プログラー角が変化した際のターゲット車両の表面電流分布をムを用いて表面電流分布の数値解析を行った。その結解析するための平面波照射とターゲット車両の条件を果、実物大のターゲット車両を対象に、準ミリ波帯レー示す。表 1 にその他の計算条件を示す。ヨー角は、ターダを照射したときの表面電流分布を実用的な計算時間

ゲット車両の中心を通り、水平面に対して垂直な直線をで解析可能であることを示すことができた。

回転軸として回転させた角度である。ターゲット車両の今後、ターゲット車両にレーダ電波を照射したときの真後ろを 0°に設定して、左右±90°の範囲で回転さ反射箇所を把握して、反射波の発生メカニズムおよびせた。

受信電力の変動メカニズムを解明する。

ターゲット車両の表面電流分布の数値解析結果を図

参考文献

3 に示す。濃い赤色の所のほど表面電流が大きく流れ[1] Eugene F. Knott, 他 , Radar Cross Section (2ndていることを示している。このように、TSUBAME3.0 とEdition), scitech publishing, (2004)

我々が開発した高精度差分スキームの FDTD（2,4）法[2] Emna Bel Kamel,

プログラムで、ターゲット車両の表面電流分布を数値解

measurements for automotive radar applications in

析することができた。
6.

他 , RCS modeling andW band, HAL, (2018)

まとめと今後の課題[3] 宇野 亨, FDTD 法による電磁界およびアンテナ

ターゲット車両からのレーダ反射波の発生メカニズム、

解析, コロナ社, p.22(2009)

および受信電力の変動メカニズムの解明に向けて、ターゲット車両の表面電流分布の数値解析に取り組んだ。

3

（様式 20）成果報告書

電流[A]

電流[A]

x

x
z

y

z

（a）ヨー角：0°y（b）ヨー角：10°y（d）ヨー角：90°

電流[A]

x
z

x
y

z

（c）ヨー角：20°

図 3. ターゲット車両の表面電流分布

4

電流[A]

（様式 20）成果報告書

TSUBAME 共同利用 令和 3 年度 産業利用 成果報告書

利用課題名 分子標的型農薬開発のためのタンパク質-阻害剤複合体シミュレーション
英文：Simulation of protein-inhibitor complex for molecularly-targeted pesticide development
利用課題責任者
田中良樹
所属
株式会社アグロデザイン・スタジオ
https://www.agrodesign.co.jp/
邦文抄録（300 字程度）
当社では医薬品で用いられている技術を応用し、分子標的農薬の研究開発を行っている。分子標的農薬の薬剤開発には、候補化合物との複合体立体構造情報が必要であるが、対象生物が多岐にわたる農薬開発においては医薬品以上に困難がある。本研究課題では、TSUBAME を利用した具体的な分子標的農薬の研究開発に繋がる情報の取得を目的とし、標的タンパク質-阻害剤複合体のフラグメント分子軌道法の条件検討を行った。

英文抄録（100 words 程度）
We are conducting research and development of molecularly targeted agrochemicals by applying thesame technology used for pharmaceuticals. The development of molecularly targeted agrochemicals
requires information on the complex structure of candidate compounds, which is more difficult than forpharmaceuticals in the development of agrochemicals that target a wide range of organisms. In thisresearch project, we investigated the conditions for fragment molecular orbital method of target
protein-inhibitor complexes to obtain information that will lead to the research and development ofspecific molecularly targeted pesticides using TSUBAME.

Keywords: フラグメント分子軌道法、農薬開発

精度の検証および PIEDA による詳細な相互作用解析

背景と目的
近年の環境保全意識の高まりにより、従来使用さを行う。

れてきた農薬の使用禁止が相次いでおり、次世代型の農薬の開発が期待されている。そこで当社では医薬品

概要

で用いられている技術を応用し、分子標的農薬の研究

次世代の農薬には、高い安全性と薬剤抵抗性への開発を行っている。分子標的農薬の薬剤開発を進める

対処が求められている。そのための開発手法の一つとには、標的タンパク質と候補化合物との複合体の立体して、特定のターゲットタンパク質を低分子農薬の標的構造を決定する必要がある。しかし、医薬品と異なり、とした分子標的農薬の研究開発を進めている。本研究

対象生物が多岐にわたる農薬開発において、その全て課題では、標的タンパク質-阻害剤複合体の分子動力の構造を実験的に明らかにするのは実現困難である。

学シミュレーションやフラグメント分子軌道法の条件検

そこで、分子シミュレーションを用いて不足している構討を行い、TSUBAME を利用した具体的な分子標的

造情報を補い、その構造に対しドッキングや分子動力農薬の研究開発に繋げる。

学、フラグメント分子軌道法といった種々の分子計算手法を適用して化合物を探索することで、化合物の改良結果および考察を加速する必要がある。本研究課題では、標的タンパ

PDB ID: 5K6R から取得されるアセト乳酸合成酵

ク質-阻害剤複合体の分子動力学シミュレーションやフ

素 (Acetolactate Sythase: ALS) とそれを阻害する

ラグメント分子軌道法の条件検討を行い、計算速度やスルホニルアミノカルボニルトリアゾリノン系除草剤の5

（様式 20）成果報告書一つである、チエンカルバゾンメチルからなる複合体間の相互作用解析を、 補因子である チア ミンニ リン酸(Thiamine Diphosphate: ThDP)、フラビンアデニンヌ
クレオチド (Flavin Adenine dinucleotide: FAD) 存
在下で FMO 計算により試みた。計算は 2 Node、Flat
MPI (28 プロセス) で行ったが、Monomer SCC 計算が終了した とこ ろで エラー によ り計 算が中 断された
(forrtl: severe (174): SIGSEGV, segmentation
fault occurred)。

まとめ、今後の課題
複数回条件を検討しつつ FMO 計算を実行したが、エラーにより計算が中断してしまったため、目標としていた計算精度の計算や分子動力学計算との比較を行
うには至らなかった。今回の実施結果を精査し、中断
原因を解明して改めて計算を行いたいと考えている。

6

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 動的均衡にある神経回路の機械学習とその場の理論的解析
英文：Machine learning in neural networks on a dynamical balance

and its field-theoretic analysis
利用課題責任者 早川 隆
Takashi Hayakawa
所属 日本大学 医学部 生体機能医学系 薬理学分野
Department of Pharmacology, School of Medicine, Nihon University
https://nu-pharmacology.com
本課題では科研費研究課題「場の理論にもとづく動的均衡にあるレザボワ計算神経回路の設計」の一環として、動物脳にみられる動的均衡という性質を持った人工神経回路のレザボワ計算と呼ばれる学習問題における性能を、数値実験と理論解析を併用し調べた。研究代表者の開発した統計力学理論により、学習時にどのよ
うな入力に対してどのような応答・情報の読み出しがどのような確率で得られるのかを見積もることができていた。そこでこの結果をもとに従来のレザボワ計算の問題点を改善する手法を提案し、TSUBAME の計算環境を活用した大規模数値シミュレーションにて学習性能が向上することを示した。本研究はリカレント神経回路のこれまで利用されてこなかった動的な性質を機械学習に応用する試みであり、上記の結果がリカレント人工神
経回路を用いた人工知能の性能を質的に向上させる手がかりとなることが期待される。
In the present study, we investigated reservoir computing with artificial neural networks whose
excitatory and inhibitory recurrent feedbacks are dynamically balanced in a similar manner to those ofanimal brain circuits. This investigation was supported by JSPS KAKENHI Grant No. JP19K20359.
Using a statistical mechanical theory that we developed in the previous fiscal year, we probabilisticallydescribed the readout from the dynamically balanced neural networks in learning. Based on thistheoretical result, we designed a learning algorithm that potentially resolves a few bottlenecks inconventional reservoir computing. With large-scale numerical simulations of learning enabled by thecomputing environment in TSUBAME, we confirmed that the algorithm improves learning
performance. The success in taking advantage of dynamical properties of recurrent neural networks inmachine learning might lead us to obtain a clue to a better design of artificial intelligence based onrecurrent neural networks.

Keywords: recurrent neural networks, GPU computing, machine learning, reservoir computing,mean-field theory
背景と目的に記憶することができると考えられ、実応用で一定の成この TSUBAME 共同利用課題は研究代表者の科功例が知られていたが、依然として扱う際の技術的難

研費課題「場の理論にもとづく動的均衡にあるレザボ

易度が高いことが指摘されており、どのような条件であワ計算神経回路の設計」の遂行を目的としたものであれば学習がうまくいくのか、理論的な解明が求められてる。研究の背景は科研費実績報告書に記載した内容ときた。

重複するものであり、これを以下に引用する。(引用開

そのような中、ランダム結合を持つリカレント神経回始）近年、人工的な神経回路を用いた情報処理技術が路については統計力学理論を用いてそのふるまいを記

社会に大きな影響を与えている。特に、入力層から出述する方法が古くから知られていた。これらのランダム

力層への処理の流れが明確に決まっているフィードフ

結合神経回路を実問題における学習に用いた場合もォワード回路は画像認識や言語処理などで高性能を発そこそこの性能を示すことが注目され、リカレント神経

揮しており、そのふるまいの理論的な解明も徐々に進回路のふるまいを理解するためのモデルケースとしてんできている。一方、自身の出力が次の時刻での入力

さかんに研究されていた[1-3]。このようなランダム結合になるリカレント神経回路を用いると過去の情報を内部

神経回路を実問題に用いる方法は「レザボワ計算」と7

（様式 20）成果報告書呼ばれる。一方、統計力学理論が適用可能な神経回

概要路と一般のリカレント神経回路には依然として隔たりが上記の背景と目的をふまえ、令和２年度までの研究あり、その中間的なケースについても理論を構築したい

で、研究代表者は動的均衡にある神経回路に外部かという要請があった。

ら入力があった場合に、神経回路がどのように応答し、そのような中、研究代表者は神経回路の結合が完

その結果どのような出力が読み出されるかを表す、（線

全にランダムではなく構造を持っている場合についての形）応答理論を構築していた。またその結果に基づき、新しいタイプの理論を世界に先駆けて構築したところで動的均衡にある神経回路を用いてレザボワ計算の学あった[4]。特に、哺乳類脳に見られるように正の結合

習性能を向上させるアルゴリズムが設計できたため、を作るニューロン集団と負の結合を作るニューロン集団

課題学習の数値シミュレーターのプロトタイプを作成しが分離しているようなケースで[図 1(a)]、回路全体の正ていた。数値シミュレーターの実行や、線形応答理論にの入力と負の入力が動的な均衡を保つような特定の条

もとづく数値解析の実行には大きな計算量が必要とさ件を満たすと、個々のニューロン活動と回路全体の集れたため、令和３年度はこの計算を TSUBAME 上で団活動が相互作用しながらパターンを形成することが複数 CPU と GPU アクセラレーターを活用して行うこと示された[図 1(b)]。とりわけ、コヒーレント状態・カオスとした。

的振動・非自明な静止状態など完全にランダム結合の神経回路では見られないような多様な挙動を示すよう

結果および考察になり、その様子が拡張された統計力学理論を用いてTSUBAME 上で動的均衡にある神経回路の学習記述できることを示せていた[図 1(a)]。そこで、この研を数値シミュレートするプログラムコードを開発し、従究代表者自身の理論を用いてレザボワ計算の枠組み

来のランダム結合神経回路に内在する問題点であを拡張することを科研費研究課題として策定した。（引る初期値依存性や時間的な記憶減衰の問題に取り用終了）組んだ。TSUBAME の GPU アクセラレーターを活用することにより、計算性能の大幅な向上を得た。

8

（様式 20）成果報告書[3]M.Lukosevicius et al., KI-Künstliche Intelligenz, 2012[4]T.Hayakawa & T.Fukai, Phys.Rev.Res., 2020
[5]L.L.Gratiet & J.Garnier, Machine Learning, 2015

まず、数値計算と線形応答理論を組み合わせた解析によりマクロ・ミクロのリアプノフ指数の上界を見積もったところ、動的な均衡にある神経回路ではミクロ
な個別ニューロン活動のレベルではカオス的であるにも関わらずマクロな集団活動のレベルでは非カオスとみなせる場合があるという観察結果を得た。また、
この上界は同じパラメータ値を持つ神経回路でも神経結合の詳細な実現値によって変化することが示せた。そこで、事前にこのリアプノフ指数がなるべく小さくなる回路を選んで学習をさせる数値シミュレーショ
ンプログラムを開発したところ、マクロな安定性を利用して学習結果の初期状態への依存性を減らせることが確認できた。
次に入力が時間的にスパースにしか与えられず、離れた時間の間の関連を学習する課題にも取り組んだ。この場合、動的な均衡にある神経回路は、前述の集団活動のリアプノフ指数がほぼゼロになるよ
うに回路結合をあらかじめ事前学習すると性能が向上することを見出した。
最後に、動的均衡にあるランダム結合神経回路の活動がガウス過程を非線形変換したものとして表されることとガウス過程の学習理論[5]をもとに、汎化誤差を抑えるための正則化入りの確率的勾配法を設計した。これによってノイズ下で有限サンプルから統計的に学習をする場合に、学習性能が向上することが数値計算によって確かめられた。
まとめ、今後の課題以上に述べたように動的均衡により定性的に異なるタイプの集団活動と個別活動のパターンが形成され、そのことが学習に貢献することが確かめられた。
科研費研究課題の補助期間は終了したが、これらの観察事実を追加の数値シミュレーションや理論解
析により強固に基礎付けし、論文発表を行う予定である。これによって、リカレント神経回路の学習、ひいては人工知能の性能を質的に向上させる手がかりが与えられたと考える。
参考文献[1]W.Maass et al., Neural Computation, 2002
[2]H.Jaeger & H.Haas, Science, 2004

9

10

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 高性能計算向け分散メモリ・ストレージ統合システムの研究
英文：The Unified System of Distributed Memory and Storages for High Performance Computing
利用課題責任者
緑川 博子
所属
成蹊大学 理工学部
http://www.ci.seikei.ac.jp/midori/

邦文抄録（300 字程度）
本研究では，スパコンをはじめとするクラスタシステムにおいて，複数の計算ノードメモリを一つの共有アドレス空間にマップし大規模共有メモリとして提供し，各ノードコアから制限のない高速アクセスを実現するソフトウエア分散共有
メモリ mSMS を構築している．これにより，単一計算ノードで広く用いられる OpenMP, OpenACC や C 言語による大規模配列への記述を，クラスタ上の大規模共有データに対して利用可能とした．本年度は，並列処理未経験である研究者に対し mSMS を提供し，音響流体解析応用 FDTD を実装して高速処理を実現し，複数ノードにまたがる
大規模データに対する複雑な境界条件の設定などが非常に容易であるなどの成果を得た．

英文抄録（100 words 程度）
The mSMS, software-based distributed shared memory system, provides a very large shared memory,which can be accessed from any CPU core on any computing-node in a cluster system without accessiblearea limitations. It realizes a parallel programming environment for shared large data on a cluster
system which is easy to use and similar to the single-node parallel programming environment, such as
C programming with OpenMP and/or OpenACC. The mSMS is used for large-scale acoustic energy
simulations in this paper and shows its effectiveness not only in the processing performance but also
in the program development.

Keywords: software distributed shared memory; partitioned global address space; cluster computing;⚫

背景と目的
クラスタ利用による高性能計算では，現在も MPI が広計算ノードメモリの範囲を超える大域共有アドレス空間の実現と計算ノード全体に配置された大規模く用いられ，分散メモリモデルによるプログラム開発の低データへの制限のない高速アクセス⚫

生産性が未だ解決されたとは言えない．これを軽減するため，PGAS（Partitioned Global Address Space）モデル

クラスタシステムにおける並列プログラム開発の生産性向上のためのプログラミン環境，API の提供と総称される様々な言語・API が提案されてきたが[1]，大域データ配列や大域インデックスを利用可能とするも本年度は，並列処理未経験である研究者に対しのはあるものの，可能なアクセス範囲がノード隣接デー

mSMS を提供し，音響流体解析応用を実装して高速処

タ領域に限らていたり，範囲を超えた大域データアクセ

理を実現し，複数ノードにまたがる大規模データに対す

スには明示的な MPI のような局所的記述が必須であっる複雑な境界条件の設定などが非常に容易であるなど

たり，定義できる大域データサイズに制限があり大規模の成果を得た.

計算には利用できないなど，多くの不自由さが存在する
概要[2][3]．

ソフトウエア分散共有メモリシステム mSMS

本研究は，スパコンをはじめとする大規模クラスタシス
テムにおいて，以下を実現する汎用かつ大規模なソフト

mSMS (multithreaded Shared Memory System)

ウエア分散共有メモリシステムを構築し，従来の並列プ

は，クラスタ上の分散メモリを大規模共有メモリとして利ログラム開発を容易にすることをめざしている[4]．

用可能にするため，各ノードのローカルメモリアドレス空11

（様式 20）成果報告書
間の他に，どのクラスタノードからもアクセスできる大規能な拡張 C 言語 MpC，（３）pragma 文 SMint に模な共有メモリアドレス空間を構築する．図１に示すようよるプログラミングと，図３に示すように，複数のプログに，各ノードが共有メモリとして提供するメモリは，「オーラミングインターフェースを提供している．

ナーページエリア」として共有メモリアドレス空間上のそ

逐次コード

れぞれ違うアドレス領域として割り当てられる．プログラ

OpenACC

OpenMP

SMint

ムが遠隔ノードの共有データアドレス領域にアクセスした場合には，ユーザ（プログラム）には見えない形で，

GPU並列

マルチコア並列マルチノード＋GPU並列マルチコア＋GPU並列マルチノード並列

SEGV シグナルで検知し，ページ単位でローカルノードに遠隔ページをフェッチして該当アドレス領域にキャッシュする． 並列プログラムにおける実行バリア同期など

マルチコア＋ノード間並列の際に，共有データの一貫性同期が行われ，各ノード

マルチノード＋マルチコア＋GPU並列にキャッシュされたページの更新情報は，そのページのホームノードに送られてページに反映される．

図 2 インクリメンタルプログラミング
逐次コードに，OpenACC, OpenMP, SMint の並列コメントを加えるだけで自由な組み合わせでプログラムが書ける

図 1 クラスタ上で大域アドレス空間を実現する
ソフトウエア分散メモリ mSMS

図3

SMS における３つの API

C と SMS ライブラリ関数利用から，SMint API まで，好みのレベルの API を利用可能

mSMS では，従来の単一ノードにおける C 言語による並列プログラム記述を，ほぼそのままの形でクラスタ上

FDTD (2,4) 法のマルチノード並列処理の大規模共有データに対して利用できるように設計され並列処理の経験がない応用分野（音響解析）の研究ている．クラスタシステムでは， マルチノード，マルチ者と共同で TSUBAME3.0 において mSMS を用いて，音CPU コア，マルチ GPU コアなど，様々な並列方式

響エネルギーシミュレーションの並列処理を行った．

が利用可能であるが，複数の API を自由に組み合

大規模な音響解析を効率的に行うために，

わせて，mSMS では高性能計算に利用できる[4]．

FDTD(2,4) 法等の高次 FDTD 法の有効性が検討され

図２に示すように，単一計算ノードで広く用いられるマている[7]．マルチノード環境でのプログラミングのルチコア向け OpenMP や GPU 向け OpenACC など

生産性を向上させるため，ページベース分散共有メと同等のマルチノード向けの API である SMint も提供

モリシステム mSMS を利用し，図４ のスケルトン

しており，逐次コードから容易に好みの API を組み合わ

コードに示す様に，複数の計算ノード全体にある大せてインクリメンタルにプログラムを作成することもできる．

域データを対象としてグローバルビューでプログラ

応用分野の計算特性や利用可能なシステム環境，

ミングを行った．mSMS を用いたマルチノード実装でユーザの好みに応じ，（１）標準の C 言語と sms ライ

は，SIGSEGV シグナルハンドラによるリモートペー

ブラリ関数による記述，（２）大域共有配列宣言が可

ジフェッチ機能による手法の他に，図 5 のスケルト

12

（様式 20）成果報告書
ンコードに示す様なリモートデータプリロード API

遠隔データ再プリフェッチ関数 sms_overload（sms_preload_array 関数）を用いて，計算利用前

ローカルノードに実データのないアドレス範囲のに事前に遠隔ノードからデータをプリフェッチするページ（図１のオーナーページ以外の領域）は，当機能もある．これを隣接ノードとの袖領域データ交

初 ， ア ク セ ス 不 可 (NO) に 設 定 さ れ て い る の で ，換に使用し，ノード内の並列化には OpenMP を利用

preload 関数で，遠隔データをアドレス空間に配置することもできる．

する際には，そのアドレス範囲のページをリードラ

FDTD 並列処理の流れ

イト可能（RW）に設定してからデータを受け取る．

図４では，sms_startup 関数で SMS が起動され各

その後，受け取ったノードで当該ページへの書き込

ノードに同一プログラムが並列に実行される．最後みがあったかどうかの有無を記録するためキャッシの sms_shutdown でマルチノード実行が終了となる．

ュページは読み出し専用(RO)に設定し直す．本応用

各ノードでは，大域データを sms_mapalloc で確保しのように遠隔データの袖領域への書き込みがない場

た後，全ノードで実行同期 sms_sync を行った後，NT合には，バリアデータ同期（sms_barrier）の際に，

時間ステップの計算ループを実行する．計算ループ

キャッシュページは廃棄（すなわちアクセス不可 NO

内では明示的な遠隔データのフェッチの記述はないへ設定）される．しかし，ステンシル計算のように，が，SMS システムにより遠隔ノードより袖データが時間ステップ毎に，繰り返し同じデータ領域を遠隔自動的にフェッチされて計算に用いられる．FDTD でノードから繰り返し読み込む場合には，キャッシュあるため，時間ループ内には２つの処理ループが含

ページは NO->RW->RO->NO のようにアクセス設定

まれ，いずれも OpenMP を用いてマルチコアによる並（mprotect システムコール）を繰り返すことになる．

列処理が行われる．各処理ループの終わりに同じ領域に遠隔データを繰り返し読み込む場合（スsms_sync_drop 関数により，実行同期と各ノードにテンシル計算のような処理の場合）には， 最初の遠おけるキャッシュページ廃棄が行われる．

隔データプリフェッチで sms_preload を用いた後，遠隔データ一括プリフェッチ関数 sms_preload

実行同期(sms_sync)のみを行いキャッシュページを

今回の FDTD や典型的なステンシル計算のように，そのまま（RW）で保持し，2 回目以降の時間ステップ

予め次にアクセスする遠隔データの範囲がわかってでは，すでに存在するキャッシュページにそのままいる応用処理では，プログラムが遠隔データにアク上書きする sms_overload を用いることで高速化がセスしてから SEGV シグナルで検知し，ページ単位に図れる．実際のコードでは，時間ブロッキングアル

データをフェッチするのではなく，予め計算開始前

ゴリズムにより，複数時間ステップ（BT）の隣接境（データアクセス前に）に次に使う連続データ領域

界データを一度に preload，overload することで，をアドレスで指定して，指定領域を含む複数ページ

遠隔データのフェッチをさらに効率化している．

を一括でプリフェッチできる sms_preload 関数を利用し，さらに高速化できる．事前に遠隔データアク結果および考察セス範囲が既知の場合，ページ毎の SEGV ハンドラがmSMS 利用による FDTD (2,4) 処理の性能評価

省き，高速化が図れる．

ここでは，前述のような高速化を測った場合

sms_preload_array 関数は，多次元大域配列の中の(Preload API 利用，ステンシル計算部分のみ) の，矩形の部分配列を指定して preload する関数で，関マルチノード環境における並列性能評価[6]を示す．

数内部で指定部分配列に含まれるページを自動的に

性能評価は，TSUBAME3.0（Intel Xeon E5-2680 v4,

計算し一括で preload する関数である．図４の２つ

14 core, 2.4GHz × 2 / node, Intel Omni-Pathの配列アップデートの for 文の前に，図５に示すよ

100Gb/s × 4, Intel MPI 2018.1.163）と 九州大

うなデータのプリフェッチを加えると，高速化が図学の ITO Subsystem A（Intel Xeon Gold 6154, 18

れる．

core, 3.0GHz × 2 / node, InfiniBand EDR 4 ×

13

（様式 20）成果報告書
100Gb/s, MVAPICH2-X 2.2）上で実行した．OpenMP ス

間を示す．図 6 のノード１の値はベースライン，すレッド数は両アーキテクチャともに 24 とする．デ

なわち，各アーキテクチャで mSMS を使用しなかっータメッシュサイズは倍精度で，z 方向に分割してた場合のシングルノードの実行時間を示している．

並列処理する．各ノードの z サイズは 1024 である．

ITO-A および TSUBAME3.0 を用いて求めたベースラ

図 6 は，1 ノードから 32 ノードまでの並列化さ

イン実行時間に対する 2-32 ノードの実行時間の比れた FDTD(2,4) ソルバーの Weak スケーリングの率は，それぞれ 1.07～1.16 および 1.09～1.30 で実行時間を示している．凡例の’Calculation’，’ある．図に示す様に，mSMS と OpenMP で並列化した

Preload’ 及び’Barrier’ はそれぞれ，計算時間，FDTD(2,4) ソルバーはほぼ理想的なスケーリング結Preload API によるデータ交換時間，バリア同期時果を示している．

図4

mSMS における FDTD (2, 4) 並列スケルトンコード

2.5Time (s)2.0

図 5 データプリフェッチ利用時のスケルトンコード

1.5
1.0
0.5
0.0
2

図 7 バリア時間

4
8
16
Number of Nodes32

エラーバーは標準偏差0.6Time (s)

0.5
0.4
0.3
0.2
0.1

図 6 FDTD(2,4)法の SMS によるマルチノード並列処理

0.0
2

weak scaling 性能 sms_preload 利用
1024x1024x1024 格子／ノード，24 スレッド／ノード

図 8 preload 時間

14

4
8
16
Number of Nodes32

エラーバーは標準偏差

（様式 20）成果報告書
また，図 7，図 8 はそれぞれ，図６の’Barrier’ 及OpenACC)とも組み合わせることが可能である．これび’Preload’ の平均実行時間を示す (エラーバー

らを組み合わせ，スパコンクラスタにおいて高性能処は標準偏差を示す)．図の様に，TSUBAME3.0 におけ理が比較的容易に可能であることを示した．また，る実験では より新しいシステムである ITO-A と比mSMS の提供する大域共有配列に対し，データの較して，比較的長い傾向があり，それぞれノード数

所在（計算ノードの内外）を意識せずに，複雑な記述の増加に伴い向上する．また，Preload API によるが可能である点は，並列処理を専門としない応用分

データ交換時間はノード数の増加に伴い飽和する傾野のユーザにとって，大きな利点であることが明らか向にあり，ほぼ理想的な結果を示している．

になった．

さらに実応用として PML 吸収境界条件を導入し，今後，応用処理プログラムを作成するユーザが容易

ヘルムホルツ共鳴器モデルの解析を行った結果を図に mSMS を利用したプログラム開発ができるように，

9 に示す[6]．このようなシミュレーションでは通常，マニュアルやサンプルプログラムなどを整備，提供す複雑な境界条件を記述する．これまでの分散データることが重要と考える．

型の MPI などでは，記述が複雑でミスを生みやすくコストが大きかった．今回，一つの共有大域データ

参考文献
M.D. Wael, et al.: “Partitioned Global Address Space
Languages”, Journal of ACM Computing Surveys (CSUR),
Vol.47, No.62, 2015
[2] 阪口裕梧，緑川博子：”グローバルビュープログラミ
ングをサポートする PGAS 言語の記述性と性能の比較”，情処学会, HPC 研究会報告，Vol 2019-HPC-170，
No 41，pp. 1-10, 2019.7[3] Y.Sakaguchi, H.Midorikawa : “Programmability and
Performance of New Global-View Programming API for
Multi-Node and Multi-Core Processing”, IEEE Pacific
Rim Conference on Communications Computers andSignal Processing, 2019.8[4] 緑川博子, 阪口裕悟:"分散共有メモリシステム mSMS
におけるマルチノード・マルチ CPU・マルチ GPU プ
ログラミング", 情報処理学会，研究報告ハイパフォ
ー マ ン ス コ ン ピ ュ ー テ ィ ン グ (HPC), Vol.2021HPC-178,No.24,pp.1-10, 2021.3[5] Ryoya Tabata, Rei Mastuda, Toshiaki Koiwaya, Sho
iwagami, Hiroko Midorikawa, Taizo Kobayashi, Kinya
Takahashi: "Three-dimensional numerical analysis of
acoustic energy absorption and generation in an air-jet
instrument based on Howe's energy corollary", The
Journal of the Acoustical Society of America, Vol.149,
No.6, pp.4000-4012, 2021.6[6] Ryoya Tabata, Hiroko Midorikawa, Ki'nya Takahashi:"Performance Evaluation of Acoustic FDTD(2,4) Method
Using Distributed Shared Memory System mSMS", 2020
International Conference on High Performance
Computing in Asia-Pacific Region HPC Asia 2020 , pp.12 , Fukuoka, Japan, Jan. 15-17, 2020. 1[7] Y. Sendo, H.Kudo, T. Kashiwa, and T. Ohtani:
The fdtd(2,4) method for highly accu-rate acoustic
analysis in three-dimensional space. Electronics andCommunications in Japan (Part III: Fundamental
Electronic Science), 86:30 – 37, 11 2003.
[1]

に対し境界条件などが容易に記述できる点は，mSMSの大きな利点として指摘された．

図9

FDTD(2,4)法へ Berenger PML 吸収境界条件(実応

用に近い条件)を導入した，ヘルムホルツ共鳴器モデルの解析結果(PML 層数: 20, 格子数: 1024^3)

まとめ，今後の課題

2020 年度から，並列処理プログラム開発経験のない音響解析研究者と共同で TSUBAME3.0 においてmSMS を利用しソフトウエアを開発し, 2021 年度，その成果が米国音響学会の論文誌に採択された[5]．
mSMS では，ユーザの経験やシステムの環境に応じた複数のプログラミング API を提供しており，既存の CPU コア，GPU コア並列 API (OpenMP,

15

16

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 廃棄物最終処分場における間隙内流体挙動の数値解析
英文：Numerical analyses of fluid dynamics in pores of landfill waste layer
利用課題責任者 鈴木和将
Kazuyuki Suzuki所属 埼玉県環境科学国際センター
Affiliation Center for Environmental Science in Saitama
URL http://www.pref.saitama.lg.jp/cess/index.html

邦文抄録
本研究では、廃棄物最終処分場内部の水やガスの流れ問題の高品質な計算スキームの開発を目的として、数値流体解析の検討を行った。具体的には支配方程式に Navier‐Stokes 方程式を用いて SUPG/PSPG 法に基づく安定
化有限要素法により離散化した。さらに,GPGPU コンピューティングによる計算の高速化を試みた。連立一次方程式の求解には GPBi‐CG 法を導入し、マルチ GPU を利用して高速並列計算によるシミュレーションを行った。
英文抄録
The purpose of this study is to establish a numerical simulation model of fluid flow in a landfill layer
with high quality and high precision. A 3D finite element method is a powerful tool for flows having
complex geometry such as porous media in landfill. However, it leads a huge amount of computation
cost. In this study, we examined to accelerate the 3D FEM by using the Graphics Processing Unit as a
general-purpose use (GPGPU).

Keywords: landfill, Numerical simulation, FEM, Navier-Stokes equations, GPGPU

ル化した流れ場において流体シミュレーションを行い、背景と目的

近年、資源開発や地球科学の分野では、岩石の流体特性と間隙形状の関係を明らかにすることを目ような多孔質媒体における流体特性を明らかにする的とする。

ため、デジタル岩石物理学の研究が進んでいる 1)。
デジタル岩石物理学は、岩石をデジタル化し、それ

概要

に対し数値シミュレーションを適用することにより、不試料には、一般廃棄物焼却施設から採取した主均質な岩石の間隙形状を考慮して水理特性といった灰と廃棄物最終処分場においてボーリングにより採物理特性を求めることが可能となる。この手法は、複取した廃棄物試料等を用いた。

雑な廃棄物層間隙形状を有する廃棄物最終処分場

これらの試料は、マイクロフォーカス X 線 CT 装置の流体挙動の解析においても有用なツールとなるもを用いて撮影した。その後、試料の CT 画像を画像のと期待される。廃棄物層内の流体挙動は、間隙形処理ソフトウェアに読み込み、モデル間隙形状の作状に強く依存すると考えられるものの、従来の方法成を行った。さらに、作成した間隙形状のデータをでは、定量的な間隙形状分布や間隙に関する幾何

STL ファイルへ変換し、この STL ファイルを基に、自学的な位置情報を表すことができなかったため、間動メッシュ生成ソフトウェア Gmsh2)を用いて有限要隙形状と流体挙動との相互作用は十分に解明され素メッシュの生成を行った。

ていない。そこで、本研究では、マイクロフォーカス X

数値シミュレーションの支配方程式には、下記

線 CT を用いて、廃棄物材料の間隙に関する 3 次元(1),(2)式に示す非圧縮性 Navier-Stokes 方程式、

ボリュームデータを取得し、このデータから定量的な連続の式を用いた。

間隙の幾何学的な情報の抽出を試みるともに、モデ17

（様式 20）成果報告書
u i
u i
p
1   u i u j 
in , (1)
+u j
=−
+
+
t
x j
x i Re x j  x j x i 
u i
= 0 inx i

ここで、f は向かい合う面を貫通する中心軸の実際に通った経路の長さ、s は f の経路の始点と終点を結ぶ直線の長さである。

. (2)

ここで、 u は流速、 p は圧力、Re は Reynolds 数、
 は計算領域である。
上記(1),(2)の支配方程式に対して、安定化有限要素法（SUPG/PSPG 法）3,4)を適用すると以下のよ
うに弱形式が導かれる。
∫ 𝑤𝑖 (
𝛺

結果および考察上部の面より垂直に流入する一様流を与え、流れの数値シミュレーションを実施した。焼却灰の間隙内
速度ベクトルを図 1 に示す。焼却灰は、コアサンプルと比較して、粒径が小さく、密に充填されているため、

𝜕𝑢𝑖
𝜕𝑢𝑖
𝜕𝑤𝑖
+ 𝑢̅𝑗
) 𝑑𝛺 − ∫
𝜕𝑡
𝜕𝑥𝑗
𝛺 𝜕𝑥𝑖

より複雑なパスと通り水が流れている様子が確認できた。

1 𝜕𝑤𝑖 𝜕𝑢𝑖 𝜕𝑢𝑗
(
+
)
𝜕𝑥𝑖
𝛺 𝑅𝑒 𝜕𝑥𝑗 𝜕𝑥𝑗

+∫
𝑛

𝑒𝑙
+ ∑𝑒=1̅̅̅𝑘∫𝛺 (𝜏𝑢

𝜕𝑢

𝜕𝑤𝑖
𝜕𝑥𝑘

)(

𝜕𝑢𝑖
𝜕𝑡

+ 𝑢̅𝑗
𝜕𝑞

𝑛

𝜕𝑢𝑖
𝜕𝑥𝑗

+

𝜕𝑢

𝜕𝑝
𝜕𝑥𝑖

) 𝑑𝛺 = 0 ,(3)

𝜕𝑢

𝜕𝑝

𝑗

𝑖

𝑒𝑙
∫𝛺 𝑞 𝜕𝑥𝑖 𝑑𝛺 + ∑𝑒=1
∫𝛺 (𝜏 𝜕𝑥 ) ( 𝜕𝑡𝑖 + 𝑢̅𝑗 𝜕𝑥 𝑖 + 𝜕𝑥 ) 𝑑𝛺 =
𝑖

𝑖

0,

(4)

ここで、 wi , q は、それぞれ式(1),(2)に対するGalerkin 項の重み関数である。また、 u i は移流速度を表し、  は安定化パラメータを表している。これら(3),(4)式に対して Crank-Nicolson 法により時間方向
図 1 焼却灰試料における数値シミュレーション結果の離散化を行うことで、有限要素方程式が導かれる。導かれた連立一次方程式の解法には GPBi-CG

そこで、これらの試料に対して定量的な間隙の構法 5)を適用した。
本研究では、計算の高速化を達成するために、

造解析を行った。中心軸に対し直交する面の中で周GPGPU を用いた並列計算を取り入れ、NVIDIA 社

囲より面積の小さくなるところ（くびれた部分）をスロの統合開発環境である CUDA とそれらのライブラリートと定義し、スロートによって囲まれた領域をポア(cuBLAS, cuSPARSE)を用いて GPBi-CG 法を適用と定義した。焼却灰とコア試料を比較すると、空隙率した。

は、焼却灰がわずかに小さかったものの、ポア数は、

また、X 線 CT 画像データからの間隙構造解析は、焼却灰が極めて多く特徴的であった。平均ポア体積、ExFact VR 及び ExFact Analysis for Porous /平均ポア有効半径は、コア試料より小さかった。また、

Particles（日本ビジュアルサイエンス）を用いて行っ

くびれた部分のスロートは、平均面積、平均有効半た。この解析により、間隙の 3 次元幾何学的情報と径ともに焼却灰が、コア試料と比べて小さかった。

して空隙中心軸（medial axis）とその過程で得られ

図 2 に、焼却灰とコア試料の間隙のネットワーク構る 連 結 経 路 （ connecting path ） あ る い は 屈 曲 度造を細線化した結果を示す（中心軸分布）。焼却灰（tortuosity）の情報を得た。空隙中心軸の決定には、

は、コア試料と比較して多くのパス（path）が存在してburn algorithm6) が採用されている。また、屈曲度

おり、比較的小さい径のポアが連結し多数のパスをは次式で定義される。

形成している。これは、流体移動が可能となるルート

𝑇𝑜𝑟𝑡𝑢𝑜𝑠𝑖𝑡𝑦 = −がより複雑になることを示唆しており、数値シミュレ𝑓𝑠ーション結果と同様の傾向であった。

18

（様式 20）成果報告書

図 2 間隙のネットワーク構造（左図：焼却灰、右図：コアサンプル）

まとめ、今後の課題

今後、間隙内流体挙動と間隙形状との関係性について評価を行い、廃棄物最終処分場の水分移動現
象のメカニズム解明につなげていきたいと考えてい
る。
参考文献
1)
2)

3)

4)

5)

6)

辻健:貯留層マネージメントに向けたデジタル岩石の利用、石油技術協会誌,84(6) 403-410 (2019)
Geuzaine, C. and Remacle, J.F.: Gmsh: a three-dimensional
finite element mesh generator with built-in pre- and postprocessing facilities, International Journal for Numerical
Methods in Engineering, 79(11), 1309-1331(2009)
Tezduyar, T.E.: Stabilized finite element formulations forincompressible flow computations, Advanced in Applied
Mechanics, 28, 1-44(1991)
Tezduyar, T.E., Mittal, S., Ray, S.E. ＆ Shih, R.:
Incompressible flow computations with stabilized bilinear
and linear equal-order-interpolation velocity-pressure
elements, Computer Methods in Applied Mechanics andEngineering, 95, 221-242(1992)
Zhang, S.L.: GPBi-CG: Generalized product-type methods
based on Bi-CG for solving nonsymmetric linear system,
SIAM J. Sci. Comput., 18, 537-551(1997)
Lindquist, W. B., et. al: Medial axis analysis of void
structure in three-dimensional tomographic images ofporous media, Journal of Geophysical Research,
101(B4),8297-8310(1996)

19

20

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 揺動電磁気シミュレーションを用いた光放射現象に関する研究
英文：A study of light emission phenomena using fluctuational electromagnetic simulation
利用課題責任者
櫻井 篤
所属
新潟大学工学部
http://www.eng.niigata-u.ac.jp/~rad/index.html
邦文抄録（300 字程度）
本 研 究 で は ， 半 導 体 発 光 素 子 （ light emitting diode ; LED ） を 用 い た 熱 フ ォ ト ニ ク ス(Thermophotonics ; TPX)発電システムに注目した．これはエミッタに LED を用いることが特徴である．LED
では，電圧を印加して非平衡状態にすることで，プランクの法則による黒体放射とは異なったエレクトロ・ルミネ
センスと呼ばれる発光を生じる．そのため，遠方場においても高効率なエネルギー輸送ができる可能性がある．
TPX システムは，大きな電力密度を有する発電手法である．また，バイアスのかけ方を工夫することで冷却機としても応用できる．
英文抄録（100 words 程度）
In this study, we focused on a thermophotonics (TPX) power generation system using
semiconductor light emitting diodes (LEDs). This is characterized by using LEDs as the emitter. In
LEDs, by applying a voltage to bring them into an unbalanced state, light emission called
electroluminescence, which is different from blackbody radiation according to Planck's law, isgenerated. Therefore, there is a possibility that highly efficient energy transportation can be
performed even in a distant field. The TPX system is a power generation method with a large power
density. It can also be applied as a cooler by devising how to apply bias.

Keywords: thermophotonics, light emitting diodes, electroluminescence, power generation
背景と目的

問題点がある．

現在，世界規模での急速な工業化に伴い，エネルギ

そ こ で 本 研 究 で は ， 半 導 体 発 光 素 子 （ light

ー消費の増加，化石燃料の枯渇や地球温暖化が大き

emitting diode ; LED ） を 用 い た 熱 フ ォ ト ニ ク ス

な問題となっている．更に東日本大震災の発生などに(Thermophotonics ; TPX)発電システムに注目した．

より，脱原子力エネルギーの潮流がある．このような現

こ れ は NTPV 発 電 シ ス テ ム と 異 な り ， エ ミ ッ タ に状から，エネルギー利用の高効率化や自然エネルギー

LED を用いることが特徴である．LED では，電圧をによって生み出される再生可能エネルギーの利用に注印加して非平衡状態にすることで，プランクの法則によ目が集まっている．

る黒体放射とは異なったエレクトロ・ルミネセンス(EL)とこれらの問題の解決策の一つとして，近接場熱光起呼ばれる発光を生じる．そのため，遠方場においても電力発電（Near-field Thermophotovoltaic ; NTPV）

高効率なエネルギー輸送ができる可能性がある．TPX

システムがある．NTPV システムでは，太陽光だけでシステムは，大きな電力密度を有する発電手法である．

なく工場の排熱など多様な熱源が利用可能である．発

また，バイアスのかけ方を工夫することで冷却機として電原理は，太陽光発電と大きな違いはないが，波長選も応用できる．

択的エミッタによって PV セルの高感度波長域に整合した熱ふく射光を放射することで高効率な発電ができる．

結果および考察

しかし，エミッタと PV セル間の距離をナノメートルスケ

TPX システムの LED-PV セル間のふく射熱伝達をールで維持する必要があり，電力密度が小さいという計算するために，揺動電磁気学シミュレーションを用い

21

（様式 20）成果報告書
た．一般的な電磁波解析手法では，有限温度を持つ物て TPX 発電システムのエネルギー輸送メカニズムにつ体から生じる電磁波を表現できない．揺動電磁気学のいて検討した．今後は異なる半導体材料の活用も視野

理論では，物体内部でランダムに熱揺動する格子振動に入れて新たな光と物質の相互作用に起因するエネル

を揺動電流としてマクスウェル方程式の電流密度の項

ギー輸送メカニズムの解明を目指す．

に導入される．従来の電磁波解析手法に揺動電磁気学の理論を組み入れることによって，2 物体間の熱ふく射解析が可能となる．
これまでの TPX に関する研究は，発電，冷却システ
ムのいずれも近接場での運用を想定しているものがほ
とんどであり，遠方場での運用を想定したシステムの研究・開発は行われていない．TPX システムを遠方場で運用するための最低条件としてシステムの発電量がLED の駆動電力を上回る必要がある．しかし，実際にFig. 1 Schematic of the TPX system.

は非放射再結合や表面フォノンポラリトンなどの非理想特性も考慮する必要がある．
そこで本研究では，TPX 発電システムの LED と PV
セルに用いる半導体材料として，２層グラフェン（bilayer graphene; BLG）を用いる．BLG はゲート電圧と構造に応じてバンドギャップが変化する材料であり，比較的低いバンドギャップを取るためバンドギャップ波
長が赤外となりやすい．通常，グラフェンは価電子帯と伝導帯がディラック点１点で交わる特異なバンド構造を持つため，バンドギャップを形成せず発光素子としては利用できない．しかし，グラフェンが２層重なった BLGの場合においては，垂直方向のゲート電圧を印加する

Fig. 2 Spectral Energy flux between the emitter

ことによって２層間の原子の反転対称性が崩れ，バンド

and the PV cell at d = 0.1 [µm].

ギャップが形成されるため発光素子として利用することができる．
計算モデルの概略図を Fig.1 に示す．この構造で解析を行った結果が Fig.2,3 である．Fig.2 より近接場領域では単層グラフェンデバイス，２層グラフェンデバイス
共に黒体放射限界を大きく上回るエネルギー伝達が起こっているのが確認できる．一方で，Fig.3 によると遠方場領域では単層グラフェンデバイスでは黒体放射下回るが，２層グラフェンデバイスではバンドギャップ波長付近で黒体放射限界以上のエネルギー伝達が達成でFig. 3 Spectral Energy flux between the emitter

きている．

and the PV cell at d = 20 [µm].
まとめ、今後の課題

本研究では，揺動電磁気学シミュレーションを用い

22

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 液体金属流れ CFD 手法の開発及び核融合研究への応用
英文：High-order flux reconstruction method for the hyperbolic formulation of the incompressible
Navier-Stokes equations
on unstructured grids
利用課題責任者
胡 長洪
所属
九州大学応用力学研究所 新エネルギー力学部門
https://www.tj.kyushu-u.ac.jp/
邦文抄録（300 字程度）
英文抄録（100 words 程度)
A high-order Flux reconstruction implementation of the hyperbolic formulation for theincompressible Navier-Stokes equation is presented. The governing equations employ Chorin’s
classical artificial compressibility (AC) formulation cast in hyperbolic form. Instead of splitting thesecond-order conservation law into two equations, one for the solution and another for the gradient,
the Navier-Stokes equation is cast into a first-order hyperbolic system of equations. Including thegradients in the AC iterative process results in a significant improvement in accuracy for thepressure, velocity, and its gradients. Furthermore, this treatment allows for taking larger timesteps since the hyperbolic formulation eliminates the restriction due to diffusion.

Keywords: Hyperbolic method, Flux reconstruction, Unstructured grid, Incompressible NavierStokes equations, Artificial compressibility method
diffusion can be significantly alleviated.
Additionally, the hyperbolic method was found toincrease the accuracy and order of accuracy ofvariables and their gradients. This can be adesirable feature in applications that require
accurate gradients.

背景と目的
High
order
methods,
such
as
flux
reconstruction (FR) method, can achieve highorder spatial accuracy on complicated geometries
using compact stencils that only involve
immediate face neighbors. Furthermore, highorder methods, such as the flux-reconstruction
and discontinuous Galerkin methods havecompact stencils which renders them particularly
suitable for computation on modern hardware
such as general purpose graphical processing
units (GPGPUs). A challenge that arises whensolving the Navier-Stokes equations using theclassic formulation is the severe time-step
restriction in diffusion dominated problems.
Even in advection dominated problems (i.e.,high Reynolds number flows), localized highdiffusion areas, either due to a high turbulent
eddy viscosity or an artificially introduced
stabilization viscosity, can have an adverse effecton stability, especially if such regions overlappedwith highly refined mesh zones. To address thislimitation in the context of the flux
reconstruction method for incompressible flows, a
hyperbolic method [1] for solving the advectiondiffusion method is implemented in the flux
reconstruction code PyFR[2]. Using hyperbolic
diffusion, the restriction on stability from

概要
This research aims to develop next generation
CFD techniques for solving incompressible, freesurface flows. Such phenomena are important ina variety of fields and applications such as
flooding and tsunami simulations, applications in
naval and marine engineering, and nuclear
fusion applications (flow of liquid metal as
plasma-facing material). This development isbased on the high-order Flux Reconstruction
method, which allows obtaining more accurate
results while utilizing modern hardware more
efficiently when compared to conventional CFD
techniques. The Tsubame super-computer wasused to carry carry out simulations for testingand validation of GPU-accelerated, high-order
CFD
code
for
solvingfreesurface
incompressible flows. The simulations aimed totest the accuracy and computational efficiency
and scaling of the code for problems involving alarge number of degrees of freedom.

23

報告書

結果および考察
Numerical
tests
of
the
hyperbolic
incompressible
flux
reconstruction
implementation showed significant reductions in
the absolute error of the field variables and thegradient of the velocity have been demonstrated .
Additionally, it has been shown that equal ordersof accuracy can be obtained for both the field
variables and velocity gradients. Numericalresults suggest that the improvement in theorder of accuracy of the velocity gradients lead toa matching improvement in the order of accuracyof pressure. Analysis shows that the time-steprequirements are significantly relaxed whenusing the hyperbolic solver.
This leads to a considerable speed-up of
convergence especially for diffusion dominatedproblems where the parabolic restriction can bequite severe.
The strong scaling performance of thedeveloped solver has been shown to be superior
to the existing INS-FR solver due to the extra
communication required for the computation of
the viscous fluxes in the latter.

Figure 2: Strong scaling performance for thedeveloped solver on multiple GPUs in
TSUBAME
Figure 2 shows the result of a strong scaling
study carried out on TSUBAME, demonstrating
the better scaaling performance of the hyperbolic
flux reconstruction method.

まとめ、今後の課題
In this report, progress on the development of ahigh-order hyperbolic method for incompressibloe
flows for use in high-accuracy and large-scale
simulations has been summarized. Satisfactory
results of benchmark tests were obtained. Future
work includes optimization of model parameters
and reducing the memory foot-print andcomputational cost of the solver.

References[1] H. Nishikawa, Y. Liu, Hyperbolic NavierStokes Method for HighReynolds-Number
Boundary Layer Flows.
[2] Witherden FD, Farrington AM, Vincent PE.

Figure 1: Convergence of the residual of thepressure for laminar flow past a sphere at
Re=100 using the classic incompressibile NS(blue) and hyperbolic method (red)PyFR: An open source framework for solving
advection–diffusion type problems on streamingarchitectures using the flux reconstruction approach.
Computer Physics Communications. 2014 Nov
1;185(11):3028-40.

Figure 1 shows the significant reduction in thenumber of iterations required by the hyperbolic
solver to reach convergence, when compared tothe classical incompressible solver.

24

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 海洋底探査を加速する AI とシミュレーション技術の開発
英文：Development of AI and simulation technologies to accelerate ocean-bottom exploration
橋本 博公
Hirotada Hashimoto
大阪府立大学 大学院工学研究科
Osaka Prefecture University
http://kyoindb.osakafu-u.ac.jp/html/110441_ja.html
海溝型地震や熱水活動観測による資源評価など，調査船を用いた海洋底探査の重要性が増している。外乱が複雑に変化する状況下では，革新的な船舶制御が求められるため，その第一歩として深層 Q 学習にもとづく
自律操船 AI の開発を行った。また，海底での採掘・揚鉱装置の開発を加速させるために，個別要素法の GPU
シミュレーションコードを開発し，その精度検証を行った。
The importance of seabed exploration using research vessels is increasing. An innovative ship
control is required in a situation where natural disturbances complicatedly change in time. As a firststep, we developed AI for autonomous ship maneuvering based on deep Q-learning. In addition, inorder to accelerate the development of mining and mining equipment, a GPGPU DEM code forsimulating seabed environment was developed.

Keywords: Seabed exploration, AI, DEM, GPU

背景と目的

ミュレーションを開発し，スパコン上で実行可能なツー

海溝型地震や海底火山噴火の予測，熱水活動観測

ルとして確立する。

による資源評価など，調査船を用いた海洋底探査の重要性が増している。気象海象が複雑に変化する状況下

概要

において，海洋底探査に課せられた高度ミッションを達

今年度は，深層 Q 学習ベースの自動避航 AI の開発成するためには，革新的な船舶制御が求められる。本を行った。ニューラルネットワークの入力には対象とす研究では，既存の制御理論では取り扱いが困難な海る海域をメッシュで分割した各セルの数値を入力する。

洋底探査における操船問題の新たな解として，深層強セルに与える数値は周辺他船との衝突予測領域と危化学習にもとづく自律操船 AI の開発を行う。スパコン

険度，および目標となるウェイポイントにもとづき算出さを利用した学習環境を構築し，学習用パラメータの最れる。空間的な情報を保持するため，畳み込みニュー

適化と分散型の強化学習を行うことで，海洋底探査のラルネットワーク層と全結合ニューラルネットワーク層を高度化と効率化に資する船舶制御を実現する。開発し用いた。

た自律操船 AI は，模型船を用いた水槽試験による検海底資源開発のための海底環境の模擬シミュレーシ

証に加えて，実船を用いた実海域での実証実験を実施

ョンについては，個別要素法の GPGPU コードを自前することにより定量的な評価を行う。これらの研究成果で開発した。 また，重力作用下の粉体挙動の計算精をもとに，今後の海洋底探査を加速させるための基盤度について検証を行った。

的な船舶制御技術として確立することを目的とする。
また，探査・開発の対象となる海底環境は，砂，泥，結果および考察

礫などの離散体で構成されており，採掘・揚鉱装置の開発した自動操船 AI を過去に実施した実船実験の開発を加速させるためには，連続体である海水との連遭遇シナリオに対して適用したところ，Fig.1 のように，成シミュレーションが求められる。本研究では，個別要輻輳する海域においても，相手船との衝突を好ましい

素法と粒子法を組み合わせた離散対・連続体の連成シ

余裕を持って回避しつつ，目的地に向かって自動航行

25

（様式 20）成果報告書できることを確認した。

両者の比較より，時間の経過とともに排出が減っていく様子が再現されていることが分かる。上部円筒内からの排出だけでなく，受け皿内のガラスビーズの盛り上がりなども良好に再現されていることが確認できる。
ただし，流量や安息角については定性的な予測精度にとどまっており，今後の精度改善が必要である。

まとめ，今後の課題
今年度は，深層 Q 学習ベースの自動操船 AI の開発と GPU シミュレーション用の DEM コードの作成を行い，両者についてバリデーションを実施した。その結果，自動操船 AI は輻輳した海域であっても目的地へと向か
いつつ，相手船との衝突を安全に回避できることを確認した。今後は，AI の改良を進めるとともに，模型実験
Fig.1

An example of ship trajectory of collision

や実船実験を通じて，その有用性を明らかにすること

avoidance by AIが求められる。
DEM の計算コードについては，検証実験と同条件開発した DEM コードの精度検証のため，排出法とでのシミュレーションの実施により，粉体挙動としては呼ばれる実験を行った。実験結果とシミュレーションの良好な精度での再現が可能であるが，流量と安息角の比較を Fig.2 に示す。

予測については課題が残る結果となった。今後は接触力のモデルを再検討したうえで，陽的 MPS 法との連成
解析コードを開発していく予定である。

Fig.2 Visual comparison of discharge flow26

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 逐次導入法による長鎖化合物の全原子自由エネルギー計算
英文：All-atom free energy calculation of long-chain compounds by chain-increment method
利用課題責任者
Nobuyuki Matubayasi
所属 大阪大学
Affiliation Osaka Univesity
URL
邦文抄録 本課題では、本グループが開発中の逐次導入法を用いて、polyethylene(PE)/polyvinylidene
difluoride(PVDF)ブレンドおよび poly(vinyl) alcohol(PVA)/polyvinylpyrrolidone(PVP)ブレンド系の混合ギ
ブス自由エネルギーを全原子分子動力学シミュレーションを用いて計算した。原子スケールのエネルギー解析
および構造解析によってブレンドの相溶/非相溶性やその原因を明らかにした。
英文抄録 In this project, Gibbs free energies of mixing of polyethylene (PE)/poly(vinylidene fluoride)(PVDF) and polyvinyl alcohol (PVA)/poly(vinylpyrrolidone) (PVP) blends were calculated by all-atom
molecular dynamics simulations using the chain-increment method being developed by our group.
An atomic-level energies and structural analyses revealed miscibility/immiscibility of the blends andtheir causes.

Keywords: all-atom molecular dynamics simulation, free energy calculation, polymer blend, mutual
miscibility

背景と目的

ンプリングが必要となり、大きな計算コストがかかる。

ポリマーブレンドとは、異なる種類のポリマーを混ぜ

本プロジェクトでは、 この問題を解決するために、合わせてできた多成分ポリマーの総称である。複合的高分子がモノマーの繰り返し構造からなる点に注目なポリマーの性質を有するポリマーブレンドは家電製品する。高分子全体ではなく、その構成要素であるモノや生体デバイス、自動車製品など多岐にわたって活用

マーの溶媒和自由エネルギーを計算することを考えされている。ブレンドの物性は、混ぜ合わせるポリマー

るのである。高分子に属するモノマー単位の溶媒和のお互いへの溶けやすさである相溶性に依存すること

自 由 エ ネ ル ギ ー は incremental chemicalが知られている。そこで、新規ポリマーブレンド材料を

potential と呼ばれ、同じ構造の繰り返しモノマーの設計する安価な手段として、計算機シミュレーションを数が多い長鎖ポリマーほど一定になることが期待さ

用いて相溶性を評価する手法の開発が期待されている。

れる。ポリマー鎖全体の自由エネルギーはブレンド系の相溶性は、熱力学的には混合前後の自incremental chemical potential の総和となるため、

由エネルギー変化である混合ギブス自由エネルギーでごく少数の incremental chemical potential 計算決まる。本研究では、ポリマー混合系における混合ギから全体の化学ポテンシャルを見積もることができ

ブス自由エネルギー評価を用いて、ポリマー溶液系のれば、計算時間を大幅に短縮することができる。

相溶性を見積もる。分子間相互作用を考慮するために
全原子分子動力学（MD）シミュレーションを行い、各ポ

概要

リマー鎖の平衡化学ポテンシャル（自由エネルギー）計

本課題では、重合度𝑁 = 100の polyethylene (PE),算から混合ギブス自由エネルギーを見積もる。しかしな

polyvinylidene difluoride (PVDF)およびそのポリマ

がら、非常に大きな内部自由度を持ち排除体積が大き

ーブレンド 系である PE/PVDF ブレンド, polyvinyl

いポリマー鎖は、自由エネルギー計算に膨大な数のサ

alcohol (PVA), polyvinylpyrrolidone (PVP)およびそ27

（様式 20）成果報告書のポリマーブレンドである PVA/PVP ブレンドについて、ぞれ溶媒 PE 分子と溶媒 PVDF 分子の寄与を表し、黄土色

モル比をさまざまに変化させた系を用意して 全原子はそれらの和を意味する。

MD シミュレーション計算を行った。さらに、得られた平衡構造に基づく incremental chemical potential 計

DmixG/N (kcal/mol)

算を行なった。incremental chemical potential から
ブレンド状態および分離状態の高分子全体の化学ポテ
ンシャルを見積もり、混合ギブス自由エネルギーを評価することでブレンドの相溶/非相溶性を判定した。さらに、
エネルギー論的解析と分子スケールの構造解析から、ブレンドの相溶/非相溶性の原因を考察した。

0.3
0.2
0.1
0
0

結果および考察

1/6

1/3 1/2 2/3 5/6
mole fraction of PVDF

1

図 1 は PE/PVDF ブレンド中の各溶質分子(PE,PVDF)

図 2:PE/PVDF ブレンド系における、モノマーあたりの混合の incremental chemical potential およびモノマー単

ギブス自由エネルギーの PVDF モル分率依存性。

位平均相互作用エネルギーの PVDF モル分率依存性を表す。左が溶質 PE の、右が溶質 PVDF の結果であ

図 2 は図 1 から見積もった PE/PVDF ブレンド系の混り、横軸の 0 が PE 単成分系,1 が PVDF 単成分系の合ギブス自由エネルギーの PVDF モル分率依存性で結果である。溶質 PE の場合、ブレンド中の PVDF のある。ブレンド状態にある全ての PVDF モル分率につ割合に対して自由エネルギーも相互作用エネルギーもいて正となっており、PE/PVDF ブレンドは非相溶であ変わらないのに対し、PVDF はブレンド中の PVDF のることが分かる。この結果は実験結果と整合している。

割合が増えるとともにどちらも負に大きくなり、化学的に

さらに、図 1 の結果から、このブレンドが非相溶である

安定化していく。これは、PE が非極性ポリマーであるのは PE と PVDF が非極性・極性ポリマーであることにがゆえに周囲に非極性分子がいようが極性分子がいよる PVDF-PE 間の親和性が弱いことが原因であるこようがエネルギー的な違いをほとんど受けないのに対とが分かった。

し、PVDF が極性ポリマーであるがゆえに同じ極性ポリマーである PVDF との結びつきが強くなるという化学的まとめ、今後の課題
全原子 MD シミュレーションを用いたポリマー相溶/

性質の違いに由来すると考えられる。
total

PE contribution<u> (kcal/mol)

Dm

incr

(kcal/mol)

solute PE

PVDF contribution

非相溶性の判定手法である逐次導入法を PE/PVDF

solute PVDF

ブレンドのような実在するポリマーを対象として実行し、

0
−0.5
−1.0
−1.5
−2.0
−2.5
−3.0
0

実験結果と整合する結果を得ることができた。これは、逐次導入法がブレンド相溶性評価に有用であることを示唆している。計算効率を向上させるために、計算時−2.0間の更なる短縮が今後の課題として挙げられる。本課−4.0

題では、十分な鎖の長さをとって鎖長 N=100 として計

−6.0
−8.0

算を行なったが、今後は混合ギブス自由エネルギーの−10.0
0 1/6 1/3 1/2 2/3 5/6 1

0 1/6 1/3 1/2 2/3 5/6 1

鎖長依存性を調べることにより、より短い長さの高分子

mole fraction of PVDF

図 1: PE/PVDF ブ レ ン ド 中 の 各 溶 質 の incremental

で相溶/非相溶性の判定が可能かどうかを調べる予定

chemical potential𝚫𝝁𝐢𝐧𝐜𝐫 およびモノマー単位平均相互作用

である。

エネルギー〈𝒖〉のブレンド中 PVDF モル分率依存性。左が溶質 PE、右が溶質 PVDF の結果を表す。赤および青はそれ

28

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 ヒト血清アルブミン-環状ペプチド複合体の相互作用解析
英文：Interaction analysis of human serum albumin-cyclic peptide complex
利用課題責任者 和久井 直樹
First name Surname Naoki Wakui
所属 独立行政法人国立高等専門学校機構 長岡工業高等専門学校
Affiliation National Institute of Technology (KOSEN), Nagaoka College
URL http://www.nagaoka-ct.ac.jp/
邦文抄録（300 字程度）
医薬品として近年注目を集めている環状ペプチドであるが、腎排出によって急速に体外へと排出されてしまう
課題を抱えている。この課題を克服するためには輸送体であるヒト血清アルブミンとの結合が重要である。環状ペプチドとアルブミンの複合体構造の相互作用を解析することで、アルブミンとの結合に重要な相互作用を明らかにし、環状ペプチド設計に役立てることが出来る。結晶構造解析によって得られたヒト血清アルブミンとコリスチンの複合体構造に対して分子動力学シミュレーションと相互作用解析を実施した。相互作用解析の結果、コリスチンの直鎖アミノ酸構造とアルブミンの極性アミノ酸の相互作用が形成されていることが明らかとなった。
英文抄録（100 words 程度）
Cyclic peptides have attracted much attention in recent years as pharmaceuticals, but they face
the problem of being rapidly eliminated from the body by renal efflux. In order to overcome thisproblem, binding to the transporter, human serum albumin, is important. Analysis of the interaction
between cyclic peptides and albumin in complex structures will reveal the interactions that areimportant for binding to albumin and will be useful for cyclic peptide design. Molecular dynamics
simulations and interaction analysis were performed on the complex structure of human serum
albumin and colistin. Interaction analysis revealed that the interaction between the linear aminoacid structure of colistin and the polar amino acids of albumin is formed.

Keywords: Cyclic peptides, Human serum albumin, Molecular dynamics

背景と目的

概要

低分子医薬品、抗体医薬品に続く第 3 の医薬品とし

結晶構造解析によって得られた、ヒト血清アルブミンて環状ペプチド医薬品が注目を集めている。しかし、環とコリスチンの複合体構造を初期構造とし、1 µs の分子状ペプチドは腎臓でろ過され、体外へと急速に排出さ動力学シミュレーションを実施した。分子動力学シミュレれてしまうという課題がある。この課題を克服するためにーションによって得られたトラジェクトリに対して相互作は血漿中に存在する輸送体であるヒト血清アルブミンと用解析を実施し、ヒト血清アルブミンとコリスチンの間にの結合が重要になってくる。ヒト血清アルブミンと環状ペ形成されている重要な相互作用の同定を行った。

プチドの複合体構造の報告数は限られており、ヒト血清
アルブミンへの結合メカニズムを解明するにはより多く結果および考察の複合体構造に対して相互作用解析を行う必要がある。

相互作用解析の結果、全シミュレーション時間中の結晶構造解析によって新たに得られたヒト血清アルブミ

30%以上にわたって相互作用を形成していたアミノ酸はンと環状ペプチドの複合体構造に対して分子動力学シ

Chain A で LYS212、ALA213、THR236、ASP324、Chain

ミュレーションおよび相互作用解析を実施し、ヒト血清ア

B では LYS281 と LYS286 であった（図 1）。LYS212 はコルブミンとの結合に重要となる相互作用を明らかにす

リスチンの直鎖アミノ酸構造上の主鎖カルボニルとの相る。

互作用を形成していた。ALA213 はコリスチンのアルキ
ル鎖との疎水性相互作用を形成していた。ASP324 は直29

（様式 20）成果報告書
鎖アミノ酸構造の末端に位置するカルボニルとの相互

参考文献

作用を形成していた。THR236 は直鎖アミノ酸構造上のIto, Sho, et al. "Structural basis for the binding

側鎖との相互作用を形成していた。Chain B の LYS281

mechanism of human serum albumin complexedと LYS286 は共にコリスチンの環構造を形成するアミノ酸

with cyclic peptide dalbavancin." Journal ofの側鎖と相互作用を形成していた。先行研究 1 のダルバ

medicinal chemistry 63.22 (2020): 14045-14053.

バンシンで見られたような環状ペプチドの環構造とアルブミンの相互作用が見られなかった。これはダルババン
シンとの結合部位が異なる事や複合体の形成様式が異なることが大きな要因であると考えらえられる。また、コリ
スチンにはダルババンシンには見られない直鎖アミノ酸構造があり、直鎖アミノ酸構造との相互作用を確保するためにアルブミンの極性アミノ酸が集中している部位に結合したと考えられる。

図 1. ヒト血清アルブミンとコリスチンの相互作用様式

まとめ、今後の課題
結晶構造解析の結果得られたヒト血清アルブミンとコリ
スチンの複合体を初期構造とし、分子動力学シミュレー
ションと相互作用解析を行った。コリスチンの構造的特徴である直鎖アミノ酸構造とヒト血清アルブミンの相互作用が形成されていることが明らかとなった。ヒト血清アル
ブミンとコリスチンの複合体構造はヒト血清アルブミンとダルババンシンの複合体構造と大きく異なっていることから、２つの複合体の形成メカニズムについて詳細な解析が必要である。

30

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 電気コイル隙間内に流れ込む冷却液挙動に関する大規模数値解析
英文：Large-scale computing of pouring coolant on the electric coil
利用課題責任者 金田 昌之
First name Surname Masayuki Kaneda
所属 大阪府立大学
Affiliation Osaka Prefecture University
URL http://www2.me.osakafu-u.ac.jp/htlab/

邦文抄録
電気コイルを冷却するための手法として用いられている，上部からノズルで冷却液を流下する現象において，冷却液が冷却対象である電気コイル内外でどのように流動するのかを解明することを目的として，フェーズフィールド法
LBM を用いた大規模二相流数値解析を実施した．コイル形状については水平角柱群を層ごとで角度をつけて積層したものを採用することで，実機において冷却液が直接流下する領域の一部を模擬した．解析結果より，角棒間隔を広げると表層の濡れ広がりが小さくなり，内部の濡れ広がりも小さくなることが分かった．これを無次元数により整理したところ，速度境界層厚さに由来した無次元数により各層の濡れ面積が整理できることが分かった．
英文抄録
To investigate the coolant behavior poured on the electric coil, the large scale two-phase simulations
were carried out for the liquid poured onto the simplified electric coil. The accumulated horizontal rod
arrays are employed for the simplified coil model, of which layer is crossed with angle. The
computational results show that the spreading area depends on the rod gap. Since the spreading on
the surface at the top layer decreases as the gap increases, the spreading in inner layer also decreases.
This tendency can be normalized at each layer by the dimensionless number.

Keywords: Cooling of electric devices, coolant behavior, two-phase flow, phase-field, lattice
Boltzmann method

背景と目的のまま解析対象とすると不正確な解析もしくは解析の

近年の電気機器の高出力化に伴い，発熱密度の上発散が懸念された．また，複雑形状を取り扱うことがで昇とその効果的な除熱が課題となっている．たとえば自き，さらに二相流解析の可能なモデルが必要であるこ

動車用のモータでは非電気伝導性の冷却液を発熱す

とも課題であった．

る電気コイル部分に直接流下することで冷却されてい本プロジェクトでは，明らかにすべき流動場と対象構る．この場合，冷却液がまんべんなくコイル内を浸潤し造に着目することで構造を簡略化することで，そこに流て流下することが求められ，そのためにノズルの位置，下した冷却液の挙動を最新の数値解析シミュレーショ

穴径，流量，穴の数などを規定することが理想である．

ン手法で解析することで，冷却液挙動を明らかにした．

しかしながら狙い通りに流体が電気コイルの隙間を通過して流下しているかは経験則に依存しており，適切概要に冷却できているかは不明である．

課題責任者の所属する研究室では，これまで格子ボ

冷却液がどのように内部を浸潤するのかを把握する

ルツマン法（LBM）を素地とした二相流解析手法に関ことができれば機器設計のヒントとなるが，実機においする研究を進めてきた．LBM は複雑形状の流動解析て直接観察することは非常に困難であることから数値が比較的容易であり，本研究室で開発してきた手法は解析による検討が有効であるとされた．しかしながら実

これまで課題であった相体積の保存性を向上させたモ際のコイル形状は大変複雑な構造をしており，これをそデルである．この手法を今回の対象に適用することで31

（様式 20）成果報告書
解析が可能となると考え，三次元二相流 LBM を用い合も同様であった．

て，簡易化した電気コイル内に浸潤・流下していく冷却液挙動に関する大規模数値解析を実施した．電気コイ
ル構造は水平角柱群を相ごとに斜めに積層したものとした．これは実際のモーターコイル構造の一部，特に流体が直接流下する箇所を取り出し簡略化したものであり，前年度と同様の形状である．コイル構造の概念図
ならびに解析ドメインを図 1 に示す．今年度は現象の対称性に鑑み，半分の領域を解析したものも実施した．これは格子解像度の検証ならびに解析高精度化のためである．

図２：角棒間隔が異なる場合の濡れ広がり

これらの因子が濡れ広がり面積に及ぼす効果を議論するために，ノズル径で無次元化した面積を，レイノ
ルズ数の平方根と角棒間隔で整理した．結果を図３に示す．いずれの層においても反比例の関係が確認できた．これはレイノルズ数の平方根が速度境界層厚さに反比例する現象と似ており，単層の場合に得られた結果とも相応する．

図１：積層構造の概念図と解析ドメイン

結果および考察
本年度は角棒間隔ならびに液体粘度（＝温度変化に相応）を変えた場合の解析を実施し，各層の濡れ面

図３：各層の無次元濡れ面積積を評価した．各層の濡れ面積は図２に示すように分割して評価した．角棒の間隔が広くなるにつれて上層まとめ、今後の課題の表層での濡れ広がりが小さくなることが小さくなり，角積層構造への二相流解析を実施してその濡れ広が

棒長手方向への濡れが顕在化することが分かる．そのり面積を整理することで，透過性壁面への流動と濡れ隙間から液体が下の層に流下するが，上層の影響を面積を明らかにした．今後は温度場解析の実装と周辺受けるため，こちらの濡れ広がり面積も小さくなること現象の解明を行う予定である．

が分かった．これは粘度を低く（＝温度をを高く）した場

32

（様式 20）成果報告書

TSUBAME 共同利用 令和 3 年度 学術利用 成果報告書

利用課題名 ヘロナミド類の抗真菌作用メカニズムの解明のための計算分子設計技術の開発
英文：In silico molecular design for understanding antifungal mechanism of heronamides
齋藤大明 / Hiroaki Saito
北陸大学 / Hokuriku University
ヘロナミド類の抗真菌作用メカニズムの解明のために、ヘロナミドを脂質膜に挿入した分子動力学シミュレーションを行い，ヘロナミドの膜内での動的構造や相互作用特性の詳細な解析を行う。さらに、化合物の膜への会合・透過・離脱の起こりやすさを定量的に評価するための自由エネルギー計算技術を確立する．
In this study, we have carried out molecular dynamics (MD) simulations of heronamides in the lipid
bilayer to understand the structure and dynamics of compounds in the membrane.

Keywords: heronamide, lipid bilayer, molecular dynamics simulation, free energy profile
背景と目的

学シミュレーション・自由エネルギー計算）を用いてヘロ

近年，化合物やペプチドを含む生体分子の膜会合のナミド類の抗真菌作用メカニズムを解明を目的に、以下

分子メカニズムの理解には，化合物やペプチドだけの３つの研究課題を行う．1. 「化合物（ヘロナミド類）をではなく，それを取り囲む脂質膜も含めた動的構造

含んだ脂質膜の構造・相互作用解析」2. 「化合物（ヘロや相互作用特性の理解が重要であるとの研究が多数

ナミド類）の膜会合の自由エネルギー計算手法の開報告され，多くの注目を集めている．抗真菌活性を発」3. 「脂質・化合物の分子構成変化に対する分子会有するヘロナミド類は飽和型の脂質膜にはタイトに合特性変化の検証」。これにより実在系の化合物・ペプ

結合する一方で，不飽和型脂質膜には弱く結合するチドの脂質膜への膜会合過程（吸着・透過・離脱）におことが報告されている．このようなヘロナミドの脂ける「動的構造」や「自由エネルギー変化」の詳細が明質膜への結合特性の違いは，ヘロナミドと脂質膜内らかになる。

における結合構造や相互作用特性の違いによるものと考えられる．このような分子会合の膜特異性の理結果および考察

解には会合分子と脂質膜の動的構造や相互作用特性

図１に本研究で取り扱うヘロナミドの分子構造を示す。

の評価が課題となるが，添加分子を含んだ混合脂質

これらヘロナミドを DMPC を 128 個，水分子を 8192 で膜は実験観測の難しさのために，これら特性は未だ構成する脂質二重層膜に１分子だけ存在する系と、明らかではなく，分子シミュレーションによる解析

20 %濃度で存在する系を作成し、各々の系での膜内が望まれている．

構造と、膜への結合特性評価を MD シミュレーションに

本研究は，化合物やペプチドの膜への特異的会合のよって評価した。

機構解明を目的として 1. 脂質の疎水鎖の長さや不飽和度を変更した脂質膜に標的化合物を挿入した
MD 計算を行い，膜内での添加分子の動的構造や相互作用特性の詳細な解析を行う．2. 高速・高精度の自由エネルギー計算法を開発し，化合物の膜への会
合・透過・離脱の起こりやすさを定量的に評価する
図 1. 8-deoxyheronamide C, heronamide C,

計算技術を確立する．

heronamide A, heronamide B の分子構造概要
本申請研究では分子シミュレーション技術（分子動力33

（様式 20）成果報告書
MD シミュレーションは全て定温・定圧条件下(T=303K,
P=1atm)で実行した。脂質の力場には CHARMM36を用い，水のモデルは TIP3 を用いた。ヘロナミドのDMPC 膜への結合特性の評価には、アンブレラサンプ
リングによる膜厚方向に対する PMF 計算によって評価した。反応座標は z 軸（膜厚方向）とし，ヘロナミド分子の水酸基を膜内中心方向に移動させた場合と，膜外方図 3. heronamide A と C の膜厚方向に対する PMF

向に移動させた場合の PFM を DMPC の膜からの離曲線

脱エネルギー曲線と定義して計算を行った。

まとめ、今後の課題(a)

今後はヘロナミド 20%濃度における膜内凝集のメカニ(b)

ズムと抗真菌機能の解明に関するシミュレーションの実施と解析を行う。

図 2. (a) heronamide C と(b) heronamide A の膜
内分子配向
図 ２ に DMPC 膜 に (a) heronamide C と (b)heronamide A を１分子挿入した系における構造のス
ナップショットを示す。heronamide C は 2 つの水酸基を膜の極性基領域に向けて配向する結果が示された
が、heronamide A の場合は、膜圧方向に対して大きくチルトして配向する様子が示され、膜内における配向特性の違いが示された。図 3 に各々のヘロナミド分子の膜圧方向に対する PMF 計算の結果を示す。自由エ
ネルギー値の比較の結果，膜からの離脱のエネルギ
ー障壁では heronamide A < heronamide C となり、
heronamide A の弱い膜結合特性が PMF 計算より示された。これらは生物活性の評価実験と対応する結果であり、シミュレーションや PMF 計算の正当性を示す結果である。さらに本研究では、ヘロナミド C の構造異性体、コレステロールやエゴステロールを 20%濃度における MD 計算も実施した。各々の系における膜内の構造や安定性に関する解析をし、それぞれヘロナミドの分子構造の違いに起因する、膜内構造や安定性の違いが解析計算により明らかとなっている。

34

（様式 20）成果報告書

TSUBAME 共同利用

令和 3 年度

学術利用

成果報告書

LRnLA アルゴリズムを用いた物理シミュレーション
英文：Simulation of Physical Processes with LRｎLA Algorithm
利用課題名

善甫 康成
Yasunari Zempo
法政大学 情報科学部
Computer and Information Sciences
http://cis.k.hosei.ac.jp/
DiamondTorre LRnLA アルゴリズムは、クロスステンシル数値計算スキームのためのデータの局所化に優れた temporal blocking アルゴリズムである。ただ計算強度は高いが、菱形のデータの鳥圧 t 会はデー
タの coalesce および aligned データアクセスは難がある。一方 FArSh データ構造は、波面型時間ブロッ
キング手続きのデータ転送のために以前導入した実績がある。そこで FArSh のインデックス付けに単純な方法を、主データストレージ、1 つの菱形形状のデータ内、および DiamondTorre 間のデータ交換においても、同様にその手順のルートを決めるために最適な方法を導くことにより、この方法をさらに発展させていく予定である。
The DiamondTorre LRnLA algorithm is a temporal blocking algorithm with good data localization
for cross-stencil numerical schemes. The Arithmetic intensity is high, however, the traversal of adiamond shape is not convenient for coalesced and aligned data access. FArSh data structure has
been introduced previously to transfer data for wavefront-type temporal blocking procedures. Here
we develop the method further by describing a simple way for FArSh indexing, as well as for spacefilling curve in the main data storage, in one diamond, and in the exchange between DiamondTorre.

Keywords: LRnLA algorithms ・temporal blocking ・data structure
coalescing の原則を満たすデータレイアウトを探すことである。
これまでの研究でキューブステンシルに合わせる
FArSh というデータレイアウトがあることを紹介した。十字型ステンシルの数値計算スキームの局在化では、2 次元シミュレーションでは菱形のタイリ
ングが[5]、3 次元シミュレーションでは正八面体元にしたタイリングが最適である。ただ、この種の形状に沿ってデータアクセスを行うことは、かなり複雑でありアクセスの coalescing が不十分である。
我々のプロジェクトの目的は、データレイアウトに 主 眼 を 置 き 、 (1) 主 デ ー タ ス ト レ ー ジ 、(2)DiamondTorre の 中 の デ ー タ レ イ ア ウ ト 、(3)DiamondTorre の間のデータ交換のためのデータ構造を開発し、配列のインデックスの計算のオバーヘッドを減少させることである。

背景と目的
物理問題の数値シミュレーション、特に波動現象
や流体力学の分野では、矩形メッシュ上のステンシ
ルスキームが非常によく使われているが、それを使ったシミュレーションコードは、性能効率が上がらないことが多い。そこで、一般的な入れ子になった
ループではなく、wavefront 型のアルゴリズムを用いると、データの再利用による効率化が図れる。また Wavefront 型アルゴリズムを用いると計算ウィンドウを作ることも可能である [1] 。
LRnLA(Locally
Recursive
non-LocallyAsynchronous) ア ル ゴ リ ズ ム [2,3] の 中 に は 、
wavefront に似たアルゴリズムもある。Torre という LRnLA アルゴリズムは、1D1Tの wavefront 型のアルゴリズム[1]と類似している。このアルゴリ
ズムは、GPU-CPU データ交換の性能効率がり、効率が高い並列計算ができる temporal blocking 法[4]の一つである。DiamondTorre LRnLA アルゴリズ
ム[3]も wavefront 型であり、十字型のステンシルに対してより良い局所性を提供する。十字型のステ
ンシルは特に FDTD のような波動現象の数値計算スキームで頻繁に使用した効率の高い計算例もある[3]。DiamondTorre を使う上で難しいのは、デー
タアクセスの局所性、ベクトル化のためのデータ・
アライメントおよび並列アクセスのための概要
先ず 2D 計算の 1 セルの幅の十字型のステンシルを例として DiamondTorre を構築から始める。
2D1T・DiamondTorre というアルゴリズムは、
菱 形 を 底 面 と す る 𝑥 − 𝑦 − 𝑡 柱 体 で 表 す 、 2D1T
temporal blocking を実現する。3Dシミュレーショ
ンの場合には、３番目の 𝑧 軸では通常のループを用いる。
35

（様式 20）成果報告書

図 1. 2D1Tで二次元の DiamondTorre シミュレーションをあらわす。wavefront と同様に領域全体を覆う。緑色で示された DiamondTorre では並列計算を行うことができる。

図 3. 主データストレージは再帰的データ構造をしている。
ピンクは 1 つの DiamondTorre で更新されるデータである。

DiamondTorre のコードの要点は、𝑡ループである。ループの繰り返しを行っている間に菱形形状の中のセルのアプデートを実施する。ループの繰り返
しで、それぞれの菱形形状は1セルに𝑥軸の右にずれていく。シミュレーション領域の右から始めると、DiamondTorre は wavefront と同様に領域全体を覆うことができる(図 1)。
(1) 主データストレージデータの初期化および出力は、同期の瞬間に可能である(図 2(1))。 DiamondTorre の菱形の底面はこのデータからロードする。なので、主データストレ
ージは菱形のタイルで構成される。𝑑次元の一般的な場合、多次元配列を使用して、配列𝐣=(𝑗1 , 𝑗2 , . . . 𝑗𝑑 ) で の 要 素 の セ ル の 座 標 は 𝐫 = 𝑗1 𝐞1 +𝑗2 𝐞2 +. . . +𝑗𝑑 𝐞d である。ここで、
𝐞1 = (3 − d, −1, −1, … , −1)
𝐞2 = (1, 1, 0, . . . , 0)𝐞3 = (1, 0, 1, . . . , 0)⋮𝐞d = (1, 0, 0, . . . , 1)である。
𝑑次元の配列は𝑍階数(モートン符号)により格納する。階数𝑅の配列は2𝑅𝑑 の要素をもつ。この配列をそれぞれの軸には細分化すると、2𝑑 の配列を取得する、その配列は2(𝑅−1)𝑑 の要素をもつ。細分化を進めると、最後には1セルになる。つまり、再帰的なデータ構造となっている(図３)。

(2) DiamondTorre の中のデータレイアウト
DiamondTorre の中のデータは GPU レジスターにロードされて、そのデータに計算を行うので、菱形の中のデータ交換ができるだけ早い方が良い。
(図 2(2))。
𝑀 × 𝑀形状の菱形の場合、1 つの CUDA スレッドに𝑀セルの計算を行うことができる(図２)。それを
𝐜 = (1, 1, . . . , 1) とする。番号が𝐣 = (𝑗1 , 𝑗2 , . . . , 𝑗𝑑 ) のセルから始まって、𝐣 + 𝑖𝐜, 𝑖 ∈ 𝑍 というルートを考えると、このルート上のセルは、𝐣B = 𝐣 − 𝑖𝑚 𝐜 のセルから 𝐣𝐹 = 𝐣 + (𝑀 − 𝑖𝑀 − 1)𝐜 のセルまで𝑀 × 𝑀の
菱 形 の 中 に あ る 。 こ こ で 、 𝑖𝑚 = min 𝑠(𝑖𝑠 ) , 𝑖𝑀 =max 𝑠(𝑖𝑠 ) , 𝑠 = 1. . 𝑑 である。{𝐣𝐹 }の集合または{𝐣𝐵 }の集合は菱形のグノモンとなっている(形状によって決まる数、つまり図形数:gnomon)。
DiamondTorre という𝑡ループの繰り返しには、菱形は一セルに𝑥軸の右にずれる。そこで、{𝐣𝐵 }のセルの出力と、{𝐣𝐹 }のセルの入力を行う。その他の
セルは菱形形状内で移動させる。
そのためデータの局所とアライメントのため、(𝑑 − 1)次元の𝑍階数配列を使用する。その配列の要
素は、𝐣 + 𝑖𝐜に対応するルートである。

図 4. 図 2（２）菱形内のデータ構造。青は{𝐣𝐵 }、オレンジ
は次の繰返しでロードして青のメモリ位置に書き換える。

(3) DiamondTorre の間のデータ交換には、以前開発した FArSh というデータ構造を使用する(図２,(３))。DiamondTorre には、GPU で保存 さ れ て い る 配列 か ら 菱形 の 底 面 の デ ータ と 、
FArSh からスロープのデータを読み込む。
実行後、上側の底面のデータを主データストレー
ジに、上のスロープのデータを FArSh にセーブし、以前保存したデータを書き換える。これにより計算
ウィンドウも可能になって、wavefront 上のデータ

図 2. DiamondTorre のために開発されたデータ構造。

DiamondTorre の菱形の底面はその配列の中の2𝑟𝑑の部分位置と一致する。セルのデータはAoS (arrayof structure) 方式で整理される。

36

（様式 20）成果報告書だけ保存している。
まとめと今後の課題
TSUBAME でも実行したことがあるアルゴリズ
ムである DiamondTorre [3]に最適なデータ構造を開発し、コードを実装した。
我々のプロジェクトにおいて開発したデータ構造を使用すると、菱形の底面をコードに実装すると非常に使いやすいものとなる。また、上述のコードで[3] に実装した2D1Tアルゴリズムだけではなく、𝑑次元の菱形の底面でも使えるようになった。十字型のステンシルを持つ数値スキームでは、𝑑 次元のDiamondTorre を実装できるようになった。𝑑次元の菱形のデータアクセス法を使用し、以前のコード[3]と比べ、更に高いパフォーマンスが得られるものと期待している。この結果は、他の数値計算法へも適用できるものと確信している。
参考文献[1] Wolfe, Michael. "Loops skewing: The
wavefront method revisited." International
Journal of Parallel Programming 15.4 (1986): 279293.
[2] Levchenko, V.D., Perepelkina, A.Y. Locally
Recursive Non-Locally Asynchronous Algorithms
for Stencil Computation. Lobachevskii J Math 39,
552–561(2018).
https://doi.org/10.1134/S1995080218040108[3] Zakirov, Andrey, Vadim Levchenko,
Anastasia Perepelkina, and Yasunari Zempo.
"High performance FDTD algorithm for GPGPU
supercomputers." In Journal of Physics:
Conference
Series,759(1),
012100.
IOP
Publishing, 2016.
[4] Endo, Toshio. "Applying Recursive Temporal
Blocking for Stencil Computations to Deeper
Memory Hierarchy." 2018 IEEE 7th Non-Volatile
Memory Systems and Applications Symposium(NVMSA). IEEE, 2018.
[5] Terrano, Anthony E. "Optimal tilings for
iterativepde
solvers."
Proceedings
2nd
Symposium on the Frontiers of Massively Parallel
Computation. IEEE Computer Society, 1988.

37

38

（様式 20）成果報告書

TSUBAME 共同利用 令和 3 年度 学術利用 成果報告書

利用課題名 GPU クラスタを用いたミリ波帯・テラヘルツ帯大規模広帯域電波伝搬シミュレーション
英文：Large-Scale Propagation Simulations of Ultra-Wideband Electromagnetic Fields in Millimeters and Terahertz
Wave Frequencies Using GPU Cluster
チャカロタイ ジェドヴィスノプ
Jerdvisanop Chakarothai
国立研究開発法人情報通信研究機構
National Institute of Information and Communications Technology
URL: www.nict.go.jp
ミリ波帯・テラヘルツ帯における電波は他の周波数帯と比較して、電波の直進性が大きく、見通し伝搬の形態で利用されることが多い。しかし、ひとたび電波伝搬経路のフレネルゾーン内に吸収体・散乱体・反射体などの障害物が存在すると、受信電力の予測手法の複雑さが増大し、不確かさが大きくなる。特に屋内伝搬環境は、壁、天井、机、棚や椅子等の様々な障害物が存在し、伝搬状況を正しく予測することは困難である。そこで、本研究では、オフィス等の室内環境においてミリ波帯（特に 28 GHz 帯）電波の受信電力を高精度に予測するために、大規模電磁界シミュレーション技術を開発した。実際に、大きさ 3 × 4.8 × 2.6 m3 規模の解析モデル(約 384
億セルで約 2 TB の総使用メモリ)に対して、大規模 GPU クラスタにより並列計算を行った結果、32 ノード(128GPUs)による計算時間は 1000 ステップ当たり、約 32 分で実用的な計算時間で終えることができた。通常の計算機では、解析が不可能であるモデルであるが、TSUBAME3.0 を用いて 28 GHz 帯屋内電波伝搬環境の解析が可能であることを示した。
Electromagnetic wave in millimeter-wave and terahertz bands is normally used in the line-of-sightcommunication systems due to its directivity since the size of the wall, building, objects, etc., is much larger than
the wavelength of the EM field. When there is an obstruction such as absorbers, scatterers or reflectors inside theFresnel zone of EM propagation path, the receiving power at the receiver fluctuates and changes abruptly.
Therefore, it is difficult to predict the receiving power in such cases, especially, an indoor propagation environment,
since it includes many scatterers such as chair, table, shelf, ceiling, wall, and so on. In this report, we havedeveloped large-scale EM simulation technique for prediction of receiving power at 28 GHz in an indoor
environment (office environment). As a first step, we have performed simulation of mock room having a size of
4.8 x 3 x 2.8 m3 (40.3 billion cells with 2 TB memory usage) by using 32 node (128 GPUs) of TSUBAME3.0 GPU
supercomputer. It takes approximately 32 minutes to calculate 1000 time steps using TSUBAME3.0. Therefore, it
has been shown that large-scale simulation of 28 GHz EM propagation inside an indoor environment is possibleusing TSUBAME3.0.

Keywords: ５つ程度環境(特に室内環境)におけるミリ波帯電波の受信電力

背景と目的
5G 無線通信や Beyond 5G の次世代無線通信技術

を高精度に予測するために、超大規模電磁界シミュレーション技術を新たに開発する。

は主にミリ波帯電波が利用される。しかしながら、ミリ波
帯電波の伝搬経路のフレネルゾーン内に吸収体・散乱

概要

体・反射体などの障害物が存在すると、受信電力の予測手法の複雑さが増大し、不確かさが大きくなる。また、本研究では、オフィス等の室内環境においてミリ波帯

波長が短いため、物体の表面形状等による散乱が生じ、（特に 5G 無線通信に使われる 28 GHz 帯）電波の受信

従来のレイトレーシング法や光学近似法による高精度

電力を高精度に予測するために、TSUBAME3.0 GPU

な解析が困難である。特にオフィス等の屋内電波環境

クラスタを用いた大規模並列化電磁界シミュレーション

においては、様々な障害物によって電波が散乱され、技術を開発した。解析に用いた計算手法は、電磁界解

複雑な電波環境を形成するため、受信電力を正確に予析でよく使われている時間領域有限差分測することが難しい。そこで、本研究では、様々な伝搬(Finite-difference time-domain, FDTD)法である[1]。

39（様式 20）成果報告書
本研究に用いたプログラムは文献[2]-[7]にも用いら
施⼯ボード

れた人体ばく露評価のためのプログラムを改良・拡張し上ランナーたものである。主な改良及び拡張した点は、以下の 3
点である。

スタッド

(1) 電気定数をセルごとに設定するのではなく、セル

2.3 mの各辺に細かく設定することができるように変更した。

zより細かく、かつ忠実にモデル構築が可能となった。

yx(2) 従来では、モデリングを行うプログラムと解析プ

下ランナー

2.4 m
3.0 m

ログラムを分離し、モデル全体を一旦ハードディスクに
図 1 模擬部屋 (⼤きさ 3 × 2.4 × 2.3 m3)

保存してから、解析の際に読み込む手順を踏んでいたが、大規模モデルの場合、ハードディスク容量が多く要したため、モデルデータをハードディスクに保存することなく、解析プログラムの中に直接組み込んだ。
(3) 様々なアンテナ(特に指向性アンテナ)を精度よく解析領域に組み込むために、アンテナ開口面の電磁界分布を波源とした解析手法に拡張した。
本研究の解析対象としては、主にオフィス等の屋内電波環境であるが、実測結果と比較検討を行うために、(a) xy ⾯

まず、模擬部屋内での電波環境を可視化することを目的とした。実際に、大きさ 3 × 4 × 2.6 m3 規模の解析モ
デル(約 384 億セル、約 2 TB の総使用メモリ)に対して、
大規模 GPU クラスタにより高速化を行った結果、32 ノード(128 GPUs)による計算時間は 1000 ステップ当たり、
約 32 分で終えることができ、実用的な時間で解析が可能であることを確認した。

結果および考察
まず、図 1 に示すように、実際に測定に用いた模擬(b) xz ⾯

部屋の大きさ 4 × 3 × 2.4 m3 の約 3/5 の体積を有する 3

図 2 模擬部屋内の電界分布

× 2.4 × 2.3 m3 の模擬部屋内に、波源としての長さ 5mm のダイポールアンテナを配置して解析を行った。模

ップ数及び計算時間ステップ間隔は、それぞれ 30,000

擬部屋は、r = 2.4,  = 0.001 S/m の施工ボード、金属

ステップ及び 1.907 ps である。並列計算に使用したノー製のランナーとスタッドから構成されており、各部の詳

ド数は 16 ノード(64 GPUs)で、使用メモリは約 757 GB で細な寸法は図 1 に示す。床面は、r = 13,  = 0.01 S/mある。30,000 ステップにかかる計算時間は、約 10 時間の厚さ 10 cm のコンクリートである。解析モデルの解像27 分で、1 ステップ当たりの計算時間は、1.255 秒であ

度は、1 mm で、8 層の完全整合層 (Perfectly matched

った。

layers)を含めた解析領域は、3,024 × 2,424 × 2,424 セ

模擬部屋内の 28 GHz における電界分布を図 2 に示ルで、総セル数は、約 164 億セルである。ダイポールア

す。図からわかるように、電界強度は自由空間に配置

ンテナの中心座標は、(x, y, z) = (2.0 m, 1.1 m, 1.3 m)でされたダイポールアンテナと同じように、アンテナ中心

あ る 。 励 振 周 波 数 は 、 28 GHz と し た 。 計 算 ス テで強く、距離が大きくなるにつれて小さくなる傾向がみ

40

（様式 20）成果報告書

扉
2.4 m

z

y
x
3.0 m
4m

図 4 模擬部屋内の電界分布
3

図 3 模擬部屋 (⼤きさ 3×4×2.4 m )てとれる。しかしながら、施工ボートやスタッド等から多重散乱波によって、単純に距離に依存して減衰するのではなく、定在波のような強弱が生じることがわかる。
また、ダイポールアンテナの周囲方向にも同様に電界強度の強弱が生じていることが分かる。
次に、図 3 に示すように、実際の模擬部屋と同じ大きさ(4 × 3 × 2.4 m3)の解析モデルを作成した。模擬部屋には、大きさ 10 × 25 × 10 mm3 の振止めが各スタッド中
心を通しており、実際の模擬部屋を完全に再現した。一図 5 ホーンアンテナを波源とした電界分布(垂直⼊射)

部施工ボードが切り取られた出入り口を設けた。扉の材質は、木材とし、比誘電率と導電率は、それぞれ 2.0
及び 0.001 S/m とした。解析領域は、3,024×4,824×
2,624 セルで、総セル数は、約 384 億セルである。並列計算に使用したノード数は 32 ノード(128 GPUs)で、使用メモリは約 1.62 TB である。8,000 ステップにかかる計算
時間は、約 3 時間 9 分で、1 ステップ当たりの計算時間は、1.422 秒であった。図 1 の解析モデルに比べて、約13%増加した。これはノード間通信が増加したことによるものであるが、実際の模擬部屋を実時間で解析することが可能であることを確認した。
図 4 に座標(x, y, z) = (2.0 m, 2.3 m, 1.4 m)に配置され

図 6 ホーンアンテナを波源とした電界分布(斜め⼊射)たダイポールアンテナによる電界分布を示す。図 4 に示る。そこで、本研究では、アンテナの開口面分布の電界すように、28 GHz における電波はほとんど施工ボード

及び磁界を波源として入力できるようにした。解析モデを透過し、スタッドによって電界レベルが低い領域(影)ルは、3 × 2.4 × 2.6 m3 とした。ホーンアンテナの開口面ができていることを確認した。また、スタッドからの散乱分布をモーメント法により求めておき、それらを波源とし波と入射波が作る定在波が確認できる。図 2 と図 4 をた。ホーンアンテナを用いた場合、開口面分布は磁流比較すると、領域が狭く、かつ開放部（扉）がないモデのみで表現できるが、任意波源の場合は電流及び磁ルではより複雑な電界分布を示していることが分かる。

流の両方をモデリングする必要がある。図 5 及び図 6

これまでは、波源としてダイポールアンテナを用いたは模擬部屋の外に波源（K バンドの標準ホーンアンテ

が、実際の 28 GHz 帯電波を発射するアンテナとしてホナ）を設置したときの電界分布を示す。図 5 では、模擬ーンアンテナ等のより指向性が鋭いアンテナが使われ部屋の壁面を平行にアンテナを配置し場合の解析結果

41

（様式 20）成果報告書である。スタッドの間を透過することが分かる。またホー

value problems involving Maxwell’s equations in

ンアンテナから細く、鋭いビームが形成されていること

isotropic media,” IEEE Trans. Antennas Propag., vol.

を確認できた。一方で、図 6 では、ホーンアンテナを反AP-14, no. 3, pp. 302–307, May 1966.

射板(完全導体)に向けて斜めに配置した場合の解析[2] J. Chakarothai, K. Wake, S. Watanabe, “Scalable

結果を示す。斜めに入射された電波は、スタッドによっGPU-parallelizedて進路が遮られているため、後ろに影ができており、複large-scale

数の細かいビームが形成されていることが見て取れる。

Applied

また、模擬部屋内の電波の反射や散乱の様子がよくわ

Journal, vol. 31, no. 6, pp. 661-668, June 2016.

かる。実際に 20,000 ステップにかかる計算時間は、約 7[3] J. Chakarothai, “Novel FDTD scheme for analysis

時間 53 分で、1 秒当たりの計算時間は、1.412 秒であっof frequency-dependent medium using fast inverse

た。従って、開口面分布を波源とした場合、ダイポール

Laplace transform and Prony’s method,” IEEE Trans.

アンテナを波源とした場合と比較して同程度の計算時

Antennas Propagat., vol. 67, no. 9, pp. 6076-6089,

間で解析が可能であることを確認した。一方、使用メモ

Sep. 2019.

リは、884 GB であった。このように、通常の計算ワーク[4] J. Chakarothai, S. Watanabe, K. Wake, “Numerical

ステーションでは、解析が不可能なものに対して、

dosimetry of electromagnetic pulse exposures using

TSUBAME3.0 の大規模 GPU クラスタを用いることによ

FDTD method,” IEEE Trans. Antennas Propagat., vol.

り、実用的な計算時間で大規模電磁界解析が可能で

66, no. 10, pp. 5397-5408, Oct. 2018.

あることを示した。

[5] チャカロタイジェドヴィスノプ, 和氣加奈⼦,

FDTD

methodelectromagnetic
Computational

fordosimetry

analysisofproblems,”

Electromagnetics(ACES)

渡辺聡⼀, 陳強, 澤⾕邦男, “超広帯域電磁界解析
まとめ、今後の課題のための周波数依存性 FDTD 法,” 電⼦情報通信

本研究では、28 GHz における室内電波伝搬環境を学会 和⽂論⽂ C, vol. J102-C, no. 5, pp. 102-113,

高精度に解析できる時間領域有限差分法を開発した。

May 2019.

本手法では、これまでのプログラムに対して、(1) セル[6] J. Chakarothai, K. Fujii, “A unified approach forの各辺に対して媒質定数を設定できるようにした。(2)treatment of frequency-dependent materials in FDTD

モデル作成のために一時的にハードディスクに解析モ

method,” ISAP 2019, Xi’an, China, Oct. 30, 2019.

デルを保存する作業を削減した。そして(3) アンテナの[7]

開口電磁界分布を波源とした解析が可能となった。例“Convergence of a single-frequency FDTD solution inとして、ダイポールアンテナ及びホーンアンテナを波源

numerical dosimetry,” IEEE Trans. Microw. Theoryとした場合の 3 × 4.8 × 2.6 m3 及び 3 × 2.4 × 2.6 m3 規模

Techn., vol. 64, no. 3, pp. 707-714, March 2016.

の模擬部屋内の 28 GHz における大規模電波伝搬の様子をシミュレーションによって可視化することができた。
今後は、解析結果と測定結果との比較を行い、解析手法の妥当性を確認し、さらに、より現実に近いオフィ
ス環境内の電波伝搬を解析し、様々な反射体・散乱体(電波散乱壁やリフレクトアレー等)を組み込んだ大規模電波伝搬解析へ拡張する予定である。

参考文献[1] K. S. Yee, “Numerical solution of initial boundary

42

J.

Chakarothai,

K.

Wake,

S.

Watanabe,（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

フェーズフィールド法に立脚した多結晶組織形成の大規模シミュレーション
Large-scale phase-field simulation on the formation of polycrystalline microstructures
三好 英輔
Eisuke Miyoshi
東京農工大学大学院工学研究院
Institute of Engineering, Tokyo University of Agriculture and Technology
https://web.tuat.ac.jp/~miyoshi/
少数の結晶粒が優先成長する異常粒成長現象の適切な予測は，高性能材料の開発において極めて重要な役割を担う．本研究では，フェーズフィールド法と複数 GPU 並列計算に基づき，異常粒成長の超大規模シミュレーション手法を開発した．さらに，数十万の結晶粒を用いた世界最大の 2D および 3D シミュレーションを行うこと
で，異常粒成長の代表的理論である平均場理論の妥当性の解明を試みた．特定の結晶粒と周囲粒との初期サイズ比，粒界エネルギー比，移動度比を変化させた系統的な大規模シミュレーションを実施し，シミュレーショ
ン結果と平均場理論による解析解とを詳細に比較した結果，特定粒の異常成長挙動（異常成長の発生の有無，
異常成長した際に到達する限界サイズ）が理論と良好に一致することが初めて定量的に明らかとなった．
This study developed a large-scale simulation method for abnormal grain growth phenomena basedon the phase-field modelling and parallel computing on multiple GPUs. The developed method wasapplied to 2D and 3D abnormal grain growth to elucidate the validity of the well-known mean-field
theory proposed by Humphreys. A series of large-scale simulations were performed while varying
the initial size ratio, boundary energy ratio, and mobility ratio of a specific grain and matrix grains.
Detailed comparisons of the simulated results and theoretical predictions demonstrated that theabnormal grain growth behaviors (i.e., whether or not the abnormal growth occurs and the limiting
size that can be reached by an abnormally growing grain) is well described by the mean-field theory.

Keywords: Polycrystals; Microstructure; Grain boundary properties; Phase-field method; GPGPU
背景と目的手法を開発し，数十万の結晶粒を用いた世界最大の多結晶固体の焼鈍過程で少数の結晶粒が優先成長

2D および 3D シミュレーションを行うことで，平均場理する異常粒成長は，集合組織や単結晶の創生におい論の妥当性を初めて明らかとすることを目的とした．

て必須の役割を担うことに加え，再結晶核生成などの
冶金学現象の発生機構としても知られ，工学的に極め概要て重要な現象である．異常粒成長を解析的に予測する

特定の結晶粒とマトリクス粒の初期サイズ比，粒界有力なモデルとして，Humphreys [Acta. Mater. 45

エネルギー比，移動度比を変化させて系統的な大規模(1997) 4231)] の平均場理論が多用されている．平均

2D・3D 粒成長シミュレーションを実施し，結果と平均場場理論では，特定の結晶粒の異常成長挙動を，その結理論による解析解とを詳細に比較した．使用した計算

晶粒と周囲のマトリクス粒とのサイズ比，粒界エネルギ

条件や平均場理論の詳細を以下に示す．

ー比，移動度比という 3 つのパラメータのみの関数とし計算条件て記述することができる．しかしながら，この理論の妥数値モデルとして，多結晶系の粒界移動を精度良く当性の詳細な検証は未だされていない．これは，理想表現可能な Steinbach ら [Physica D 134 (1999)

的なモデル系に対する実験の困難さと，従来の数値研385] のマルチフェーズフィールド (MPF) モデルを用

究における計算スケールおよび計算精度の限界に強くいた．計算領域サイズは，2D で 12,2882x2（x = 1:

起因するものである．

差分格子間隔），3D で 12803x3 とし，時間増分 Dt は，本研究では，フェーズフィールド法と複数 GPU 並列

陽解法の安定条件を満たす範囲でなるべく大きく設定

計算に基づき，異常粒成長の超大規模シミュレーション

43

（様式 20）成果報告書

Fig. 1 3D polycrystalline system for simulating the abnormal growth of a specific grain: (a)initial state, (b) 100,000 steps. Different grains are distinguished by different colors. The matrix
grains are visualized only for a half of the entire domain to make the specific grain visible.
した． Figure 1 (a)に例示するように，粒径 R，粒界エネルギー  ，モビリティ M の球形の結晶粒（特定粒，specific grain）を領域中央に配置し，その周囲を平均
粒径<R>，粒界エネルギー<>，モビリティ<M>の約
30 万の結晶粒群（マトリクス）で充填した．2D 計算には16GPUs，3D 計算には 64GPUs を使用した．
以上の条件の下，特定粒とマトリクスとの初期粒経
比ini = Rini /<R>ini，粒界エネルギー比 =  /<>，モ
ビリティ比 = M /<M>を変化させて計算を行い，特定粒の異常成長挙動を下記の平均場理論と比較した．
平均場理論

Fig. 2 Abnormal growth conditions of the specific

Humphreys は，Hillert の平均場解析に立脚し，マgrain shown in Fig. 1 as a function of the grain

トリクス中の特定粒の異常成長挙動を，粒経比 = R /size ratio  = R / <R>, boundary energy ratio  = <R>，粒界エネルギー比 =  /<>，モビリティ比 = M/<>, and mobility ratio  = M /<M>.

/<M>の関数として記述する解析的モデルを提示した．
この理論によれば，特定粒の異常成長が生じる最小粒

る．が実線の上側（max）に達したところで異常成長は

経比min，異常成長した粒が到達する最大粒径比 max

止まり，以後， は一定値max を保つと予想できる．

は次のように表される．

min = 2 − 2  (  −  ) ,結果および考察max = 2 + 2  (  −  ) .

Figure 1(b)は，シミュレートされた異常粒成長挙動の一例として，粒界エネルギー比 = ，モビリティ比

この関係を図示すると，Fig. 2 のようになる．特定粒と = 4 の条件における 100,000 steps 時の組織形態のマトリクスとの粒経比 が実線の下側（min）より上にあを示す．特定粒が周囲のマトリクスを蚕食して優先成長

れば，特定粒は異常成長し，粒径比 が次第に増加す44

（様式 20）成果報告書

Fig. 3 Temporal variations in size ratio  as functions of boundary energy ratio  and mobility
ratio  compared to the theoretical diagram shown in Fig. 2: (a) 2D and (b) 3D simulations. Plots
with different colors indicate results for different initial size ratios ini = Rini /<R>ini.

する，異常粒成長の典型的な描像が得られている．

長挙動（異常成長の発生の有無や最終到達サイズ）は，

Figure 3 は，計算と平均場理論との比較の例として，

平均場理論による予測と良好に一致することが示され

粒界エネルギー比 =  /<>を 0.75～2，モビリティ比

た．以上の内容に加え，より詳細・広範なデータや，計= M /<M>を 0.75～4 と変化させ，さまざまな初期粒径

算と理論との若干のずれの原因に関する考察を盛り込比から計算を行い，得られた粒径  の時間変化を理論

んで論文化し，現在，国際誌に投稿中である．

線図（Fig. 2）上にプロットした結果を示す．ここで，異なる色のプロットにより，異なる初期粒径比 ini = Rini

まとめ、今後の課題/<R>ini からスタートした計算結果を区別している．また，フェーズフィールド法と複数 GPU 並列計算を活用しプロットの近くに示した矢印は，粒径比  の時間変化のた大規模シミュレーションにより，Humphreys の著名

方向（上矢印：増大 [d/dt > 0]，下矢印：減少 [d/dt

な異常粒成長理論の妥当性を初めて定量的に明らか< 0]）を示す．Figure 3 より，異常粒成長の発生の有無とした．また，この成果を通じて，異常粒成長等の多結（d/dt の正負）は，ほとんどの結果において平均場理

晶組織形成における粒界物性（エネルギー，モビリテ論による予測と一致することがわかる．さらに，異常成ィ）の重要性が改めて浮き彫りとなった．今後は，分子長発生時の最終的な粒径比は，同理論の解析解の近動力学による原子論的解析と，データ同化をはじめと傍に収束している．

するデータ科学手法を援用することで，情報が著しく不これらの結果から，マトリクス中の特定粒の異常成足している状況にある各種粒界物性値の評価・蓄積の45

（様式 20）成果報告書
方法論を確立し，大規模フェーズフィールド計算に立脚した多結晶組織形成予測の更なる高度化を目指す．

46

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

水圏動物が有する酵素の構造解析
Structural analysis of enzymes from aquatic animals
潮 秀樹
Hideki Ushio
東京大学大学院農学生命科学研究科
Graduate School of Agricultural and Life Sciences, The University of Tokyo
http://www.a.u-tokyo.ac.jp/index.html
長鎖高価不飽和脂肪酸の生合成能が高いフトメリタヨコエビ（Melita rylovae）に着目し，トランスクリプトーム
解析を行なった．トランスクリプトームから脂肪酸伸長酵素 Elovl6 を同定し，演繹的にアミノ酸配列を得た．得られたアミノ酸配列をもとに立体構造予測を行った．予測された立体構造について Ramachandran plot の描画などによって検証した結果，活性部位を含む大部分は信頼度が高いことが示された．これを用いて分子動力学
シミュレーションを行い，M. rylovae の Elovl6 の基質特異性や温度依存性について検討している．
Transcriptomic analysis of Melita rylovae identified the expression of Elovl6. The amino acidsequence of M. rylovae Elovl6 was determined by translation of the identified mRNA and used to
predict its structure. The predicted structure was evaluated by the Ramachandran plot. The results
showed that most of the predicted structure of the M. rylovae Elovl6 including active site has highconfidence. Molecular dynamics simulation by using the predicted structure has been running toinvestigate the specificity or temperature-dependent activity of the enzyme.

Keywords: aquatic bioscience, LC-PUFA, Elovl, amphipoda, crustacea

背景と目的

行った．

水圏動物の組織で発現する酵素には，水圏動物に特有なものや陸上動物よりも活性が高いものなどがあ結果および考察る．これらの酵素の反応機序を明らかにすることによっヨコエビ類の数種について LC-PUFA 生合成能を比て，遺伝子工学的に活性の向上や基質特異性の改変

較した結果，フトメリタヨコエビ（Melita rylovae）が高いが可能となり，様々な産業利用が拓かれる期待される．

LC-PUFA 生合成能を有することが示唆された．そこで，

しかし，それらの詳細な反応機構については未解明な

M. rylovae の Elovl6 について，立体構造予測を行な部分が多く，ゲノム編集による活性の向上や基質特異った（図 1）． Ramachandran plot の描画，及び立体性の改変において大きな障壁となっている．本課題で構造既知のヒト由来 ELOVL7 の立体構造との比較には，トランスクリプトーム解析で得た mRNA 配列から演

よって，予測構造の妥当性を検証した．その結果，概ね繹的に推定したアミノ酸配列をもとに立体構造の予測信頼度が高い立体構造が得られた．ただし，C 末端，Nおよび分子動力学シミュレーションを行い，水圏動物由末端の予測構造については信頼度が低かった．得られ来酵素の反応機構について明らかにする．

た立体構造をもとに，力場として charmm36m を用い，
Gromacs による分子動力学シミュレーションを進めて

概要

いる．

本課題ではヨコエビ類を対象とし，長鎖高価不飽和
脂肪酸（long-chain polyunsaturated fatty acids,
LC-PUFA）伸長酵素（elongation of very-long-chain
fatty acids, Elovl）の発現についてトランスクリプトー
ム解析によって調べ，アミノ酸配列を演繹的に推定した．
得られた Elovl のアミノ酸配列から立体構造の予測を

図 1 Elovl6 の予測構造

47

（様式 20）成果報告書
まとめ、今後の課題予測された Elovl6 の立体構造は一部で信頼度が低かったが，活性部位を含む大部分では高い信頼度であった．今後は，M. rylovae の Elovl6 の基質特異性や温度依存性について検討を行う．また，翻訳後修飾の有無などについても調べる必要がある．

48

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 蛋白質天然変性領域の動力学特性に関する計算科学的検討
英文：Computational study on dynamics of intrinsically disordered proteins
利用課題責任者
Kota Kasahara
所属
立命館大学生命科学部

独自の拡張アンサンブル分子動力学法である Virtual-system canonical coupled molecular dynamics 法（VcMD）の改良を行い、転写因子 TGIF-1 の天然変性領域（Intrinsically disordered region; IDR）の動的特性の解析を行った。VcMD の改良により反復計算の過程で得られたすべてのアンサンブルを reweighting
し、カノニカル分布を近似的に得ることができるようになった。TGIF-1 では IDR がリン酸化を受けることで構造ドメインの DNA 認識ヘリックス上の Arg と塩橋を形成し、DNA 結合を競争的に阻害するメカニズムが示唆された。
We have developed a new version of our original extended-ensemble molecular dynamics method,
named virtual-systfem coupled canonical molecular dynamics (VcMD). In this method, all theconformational ensembles obtained in the iterative simulations can be reweighted to estimate thecanonical ensemble. This method was applied to investigate dynamic properties of an intrinsicallydisordered region (IDR) of the transcription factor, TGIF-1. As a result, representative structures inthe simulated ensemble show salt-bridge formation between DNA recognition helix andphosphorylated IDR, suggesting that interactions between them competitively inhibit recognition of
DNA recognition.

Keywords: ５つ程度

背景と目的手法であり、本プロジェクトでは IDRs の広大な構造空間

蛋 白 質 の 天 然 変 性 領 域 （ intrinsically disordered

を効率よくサンプリングするために方法論のさらなる改regions; IDRs）は特定の立体構造を持たないにも関わら良を行った。これを応用し、転写因子 TGIF-1 の IDR がず様々な重要な機能を担っている。その構造は常に揺リン酸化を受けて DNA 結合ドメインの機能を制御するメ

れ動いているため、従来の構造生物学的手法ではそのカニズムを調べた。IDR の一部をペプチドとしてモデル

分子機構を解析することが困難である。これに対し、分化し、DNA 結合ドメインとの相互作用を解析した。IDR

子動力学法(molecular dynamics; MD)による動力学特についてはリン酸化状態と非リン酸化状態の 2 状態の性の解析が、分子機構解明に向けた有効な手段となる。

モデルについて計算を行った。

一方で、IDRs の極めて高い構造自由度は、サンプリン
グすべき構造アンサンブルを膨大なものとするため、結果および考察

MD による計算も容易ではない。そこで本プロジェクトで

VcMD 法 で は 任 意 の 構 造 パ ラ メ ー タ （ collective

は独自の手法である Virtual-system coupled canonicalvariables; CVs）を設定し、バイアスポテンシャルによっmolecular dynamics 法（VcMD）のさらなる開発を行い、て CVs に沿った構造変化を促す。バイアスポテンシャル

この手法を活用することで IDRs の動力学特性の解明をは直接的に構造変化を促すように印加されるのではな

目指す。

く、望まない構造変化を起こした際にそれを妨げるように働く。従来法のように人工的なポテンシャルによって概要

構造変化を積極的に引き起こすのでなく、望ましい構造

VcMD は我々のグループで考案し開発を行ってきた変化が起こる確率を高めるという方針で構造変化を促
49

（様式 20）成果報告書す。このため得られたアンサンブルには人工的な構造

べた。

変化（ゆがみ）が少なく、より自然に近いアンサンブルがその結果得られた構造アンサンブルより、特に安定だ得られるものと期待できる。

った代表構造を図 1 に示す。リン酸化状態では IDR が本研究ではさらなる方法論の開発を進め、状態密度

認識ヘリックスと直接相互作用している。特に負電荷を分布の推定方法などについて改良を行った。VcMD 法

帯びたリン酸期が、kink 周辺の Arg 残基と塩橋を形成

では反復計算によって状態密度分布を推定するが、反していた。この部位は本来負電荷を帯びた DNA と相互復計算の過程で得られたアンサンブルは最終的には用作用するべき部位であり、リン酸化した IDR は DNA と競いることはできず、反復計算で得られた状態密度を固定

争的に結合するものと考えられる。一方非リン酸化体でして production run を改めて実行する必要があった。こは IDR は認識ヘリックスとは逆側の表面へと接触しておの点を改良し、反復計算で得られたすべてのアンサン

り、DNA 結合に干渉しないことが示唆された。

ブルを再重み付けしてカノニカルアンサンブルを近似

すなわち TGIF-1 では IDR におけるリン酸化、脱リン的に算出する方法論を考案し、実装を行った。

酸化サイクルによって、DNA 結合ドメインと IDR の相互

次に VcMD 法を用いて転写因子 TGIF-1 における作用が制御されているものと考えられる。すなわち、IDR

IDR の動力学特性の解析を行った。TGIF-1 はホメオドの動的な構造変化によって転写因子の機能制御が行メインによって DNA を認識する転写因子であり、がんな

われていることが示唆された。

どの疾患との関連も高く医学・薬学の観点から重要な転写因子である。ホメオドメイン近傍に IDR を有しており、まとめ、今後の課題

IDR の存在によって DNA 結合能が制御されることが示

独自の拡張アンサンブル法である VcMD の改良と、転唆されている。特に IDR がリン酸化されることで転写活

写因子 TGIF-1 における IDR の動的特性解析への応性が阻害される。本研究では TGIF-1 ホメオドメインおよ用を行った。リン酸化が IDR の動的特性を変化させ、転び IDR の一部を切り出したペプチドの 2 分子を含むモ

写因子の機能制御を行うメカニズムが示唆された。一方デルを構築し、ホメオドメインとペプチドの距離を CV とで IDR の広大な構造空間のうちどの程度がサンプリング

して設定して計算を行うことで、IDR がホメオドメインに対されたのか、その解釈や妥当性の評価は不十分である。

してどのように作用するのか、その分子メカニズムを調今後さらなる詳細な解析を行っていくことが望まれる。

べた。特に IDR についてリン酸化状態と非リン酸化状態の両方について計算を行い、リン酸化による影響を調

図 1. 計算より得られた代表構造。(左)リン酸化状態、(右)非リン酸化状態。DNA 結合ドメイン（ホメオドメイン）を緑
で、IDR を青で示した。図中右側の kink したヘリックスが認識ヘリックスである。

50

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

超流動ヘリウムにおける量子乱流の数値的研究
Numerical study of quantum turbulence in superfluid helium
湯井 悟志
Satoshi Yui
慶應義塾大学 自然科学研究教育センター
Research and Education Center for Natural Sciences, Keio University
https://www.sci.keio.ac.jp/
我々は，量子流体力学における重要問題である 2 流体結合ダイナミクスを研究している．本課題では，量子乱流と常流体乱流の結合ダイナミクスを研究した．そのような 2 流体同時乱流は，未解決問題である量子乱流の T2 状態の原因だと予想されてきた．我々は，2 流体結合ダイナミクスの数値計算を行い，2 流体同時乱流を解析した．我々の研究結果は常流体乱流が量子乱流を強化することを明らかにし，得られた結果は実験と一致していた．
We study the coupled dynamics of two fluids, which is an important problem in quantum
hydrodynamics. In this project, we studied the coupled dynamics of quantum turbulence andnormal-fluid turbulence. The dual-turbulent state has been expected to be the cause of the unsolved
T2 state of quantum turbulence. We performed numerical simulations of the coupled dynamics of
two fluids, and analyzed the dual-turbulent state. The results of our study showed that thenormal-fluid turbulence enhances the quantum turbulence, and the obtained results are consistent
with the experiments.

Keywords: 量子流体力学，超流動，量子渦，量子乱流，2 流体モデル，渦糸モデル

2 流体の両方が乱流である状態は，量子乱流におけ背景と目的我々の研究目標は，超流動ヘリウム中の量子乱流の物る長年の未解決問題である T1-T2 遷移の原因であると理を明らかにすることである．液体ヘリウム 4 は，大気圧予想されてきたが[4]，2 流体結合ダイナミクスを考慮し下の温度 2.17 K 以下で粘性が消失する超流動性をなければ取り組めない問題である．ここで，T1-T2 遷移示す．ティサとランダウの 2 流体モデルによると，超流動状とは超流体が T1 状態と T2 状態の 2 段階で乱流遷移

態の液体ヘリウム 4 は非粘性の超流体と粘性をもつ常する現象のことである．これに対して，T2 状態は 2 流体

流体の混合流体として理解できる[1,2]．本課題で取り両方が乱流の状態である，という仮説がある．我々は，組む量子乱流とは超流体の乱流のことであり，量子渦

量子乱流の T2 状態を解明するために，2 流体モデルの(渦芯まわりの速度循環が量子化された渦)のタングルと数値計算を用いて 2 流体同時乱流状態の研究を行なして実現される[3]．量子乱流は，他の様々な系(超流った．具体的には，超流体の渦糸モデルと常流体のナビ

動ヘリウム 3，原子気体ボース・アインシュタイン凝縮体，

エ・ストークス方程式を連立して数値計算を実行する．

中性子星など)においても出現し得る．したがって，量子

我々の先行研究において，2 流体結合ダイナミクスの乱流の研究は，物理学の広範囲の分野に対して重要計算の基本的な導入には成功し，興味深い物理を明な意味をもつだろう．

らかにしてきた[5,6]．しかし，T2 状態の量子乱流等に本課題では，超流動ヘリウム 4 における超流体と常流

対しては数値計算コストが非常に大きい等の問題点が体の結合ダイナミクスの解明に取り組んだ．量子乱流状あった．そこで，渦糸モデルの数値計算に対して Fast

態では，量子渦を介して 2 流体間に相互摩擦力が働く

Multipole Method [7,8]などの計算高速化を適用しようになり，2 流体が互いに影響を与えながら運動を行う．

て TSUBAME の性能を活用する．

そのような 2 流体結合ダイナミクスは，量子流体力学における根本的かつ重要なテーマである．
51

（様式 20）成果報告書
概要の値が一定値の周りを揺らぐようになる．これは，エネルギー

超流動とは非粘性の流れのことであり，極低温において量の注入と散逸がつりあった統計的定常状態に到達したこと

子力学的な効果がマクロに現れた現象である．そのようなを示している．

超流動の非線形・非平衡状態である量子乱流は，超流

以上のようにして得られた量子乱流および常流体乱流

動ヘリウム 3 および 4，冷却原子気体，中性子星等分野をの統計量を解析していく．統計量は，統計的定常状態に横断する重要な現象である．我々は，典型的な系である

おいて時間平均されている．量子乱流の統計量は上述の超流動ヘリウム 4 の量子乱流の研究を行う．本課題では，

渦糸長密度 L であるが，常流体の統計量としては常流

量子乱流の長年の未解明問題である T1-T2 遷移に関連

体速度ゆらぎ I を採用する．常流体速度ゆらぎ I は，した研究を行う．

量子乱流の可視化実験等により調べられており，それらの実験と値を比較することが可能である．数値計算では，外結果および考察

力の振幅を調整することにより，常流体乱流のゆらぎの大

超流体に対する渦糸モデルでは，量子渦の芯の太さを無

きさを変化させることができる．各平均常流体流速 Vn に対視して糸で近似する[3]．通常の古典流体と違い，超流動

して，外力の振幅を変化させていくつかの数値計算を実行

ヘリウム 4 に対しては渦糸モデルが非常に有効な方法としてし，ゆらぎ I と渦糸長密度 L の関係を探索した．さらに，

知られている．超流動速度場は，全渦糸に沿って積分を

ゆらぎ I が実験と一致するように外力の振幅を調整した

行なうビオ・サバール則によって決定される．温度 0 K にお状況において，平均流速 Vn と渦糸長密度 L の関係をいては，渦糸は超流動速度場に完全に追随して運動を行解析した．その結果，多くの実験で知られているように[3,4]，

なうとされる．しかし，有限温度においては常流体との相互L1/2 が Vn に比例することを確認した．この比例関係の係摩擦の影響も受けて運動を行なう．一方，常流体はナビ数γは，量子乱流の T1 状態と T2 状態を判別するパラメ

エ・ストークス方程式のような運動方程式に従い，外力としータであり，T2 状態のγは T1 状態のものより大きい傾向て量子乱流から受ける相互摩擦力が考慮される．このようがある．我々の数値計算の結果，2 流体同時乱流状態のに 2 流体は相互摩擦を介して互いに影響を与えながら運γの値は，T2 状態のものと一致することが明らかになった．

動 を 行 な う ． 今 回 の 研 究 で は ， 常 流 体 に Arnold

以上の結果は，背景と目的で説明した「T2 状態は 2 流体

–Beltrami–Childress 型の外力をさらに与えることで常の同時乱流状態である」という考えを支持している．

流体を乱流にしている．また，常流体には非圧縮条件を課まとめ、今後の課題

している．
mm)3 の空間 3 次元の系で行

以上のように，我々は量子流体力学における長年の未なう．x, y, および z 方向に周期的境界条件が課されてい解明問題である T1-T2 遷移に関連する研究を行なった．

る．x 軸方向に熱対向流を印加して量子乱流を駆動する．

具体的には，T2 状態では常流体と超流体の両方が乱温度 1.9 K に対応するパラメータを用いる．超流体の時間流である，という仮説を確かめるために，2 流体結合ダイ

積分は 4 次の Runge–Kutta 法，常流体の時間積分はナミクスの数値計算を実行した．結果として，2 流体同

2 次の Adams–Bashforth 法で行う．平均の常流体流

時乱流の統計的定常状態を得ることに成功し，得られ速を Vn = 2.5 mm/s, 3.0 mm/s, 4.0 mm/s, 4.5 mm/s

た統計量は T2 状態の実験値と一致していた．これらのとして数値計算を行う．

結果は，T2 状態の仮説を支持しており，大きなインパク

数値計算は，体積(1.0

トがあるだろう．

数値計算の結果，2 流体同時乱流の統計的定常状態が得られた． 量子乱流の統計量としては，渦糸長密度 L

しかしながら，研究結果を論文としてまとめて投稿した(単位体積あたりに含まれる渦糸の長さ)に注目した．初期

ところ[9]，いくつかの批判を受けて掲載拒否となった．批状態として少数の渦糸を配置しているが，熱対向流からエ判の一つは，数値計算を行なったパラメータの範囲に関

ネルギー注入を受けて時間とともに L の値が増加していく．

することであった．そこで，我々は TSUBAME の計算機

今回の数値計算の場合，時間が 5 s 程度経過すると L資源を活用してさらに広いパラメータでのデータを集める

52

（様式 20）成果報告書など，研究内容の説得力を上げていたが，利用期間中には計算が完了しなかった．今後も，TSUBAME を利用してデータを集め，近い内に論文として結果を発表したいと考えている．
参考文献[1] L. Tisza, Nature 141, 913 (1938).
[2] L. Landau, J. Phys. U.S.S.R. 5, 71 (1941).
[3] R. J. Donnelly, Quantized Vortices in Helium

II, (Cambridge University Press, Cambridge,
England, 1991).
[4] J. T. Tough, Superfluid Turbulence, in Prog.
in Low Temp. Phys. 8, Chap. 3 (1982).
[5] S. Yui, M. Tsubota, and H. Kobayashi, Phys.
Rev. Lett. 120, 155301 (2018).
[6] S. Yui, H. Kobayashi, and M. Tsubota, Phys.
Rev. Lett. 124, 155301 (2020).
[7] R. Yokota, T. Sheel, and S. Obi, J. Comput.
Phys. 226, 1589 (2007).
[8] R. Yokota et al., Comput. Phys. Commun. 180,
2066 (2009).
[9] S. Yui et al., arXiv:2105:09499.

53

54

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 新規 BRCAness 誘導阻害剤開発に向けた理論研究
英文：Theoretical study for the development of novel BRCAness inhibitor
小関 準
Jun Koseki
国立大学法人東海国立大学機構名古屋大学大学院医学系研究科
Graduate School of Medicine, Nagoya University
URL
邦文抄録（300 字程度）
近年、合成致死の原理を利⽤した PARP 阻害剤による BRCA 変異陽性がんの治療法が再脚光を浴びている。
PARP 阻害剤の有効性は様々な癌腫においても確認されているため、幅広くリバース・イノベーションを起こしうる可能性がある。そこで、未だ詳細が明らかではない「BRCAness」に関わる分子を標的とした、革新な新規誘導阻害剤の研究開発を行うことを目的とし、研究を推進した。共同研究先では現在、「BRCAness」に関わる有効分子標的探索を並行して実施しているため、TSUBAME を用いた創薬標的分子解析では、まずは Alphafold2 を用いてヒトの全タンパク質の立体構造予測を実施した。
英文抄録（100 words 程度）
Recently, PARP inhibitors, which utilize the principle of synthetic lethality, are gaining renewed
attention for the treatment of BRCA mutation-positive cancers, and their efficacy has been confirmed
in various types of carcinomas, which may lead to a wide range of reverse innovation. Therefore, we
promoted research aimed at developing innovative novel inducible inhibitors targeting molecules
involved in "BRCAness," the details of which are still unclear. Since our collaborators are currently
searching for effective molecular targets related to "BRCAness" in parallel, we first used Alphafold2 to
predict the 3D structures of all human proteins in the analysis of drug target molecules.

Keywords: ５つ程度 BRCA, Drug Terget, Drug Design, Structure Prediction, Structure Analysis

背景と目的

概要

近年、合成致死の原理を利用した PARP 阻害剤によ本プロジェクトでは、Domainfocused CRISPR スク

る BRCA 変異陽性がんの治療法が再脚光を浴びておリーニングから得られた合成致死性を示す機能的ドメイ

り、国内外の臨床試験から様々ながん腫における有効ン近傍のタンパク質対する創薬に備えて、Alphafold2性が実証されてきた。2016 年に転移性前立腺がんのを用いてヒトの全タンパク質の立体構造を予測し、さら大規模解析から、標的候補として ATM を含む 15 個のに申請者が作成した独自のプログラムにより、生体内遺伝子に変異が同定されている。しかしながら BRCA

環境を模倣した分子動力学シミュレーション用のインプ以外の変異は 1%前後と非常にレアである ことや、

ットファイルを自動生成した後、エネルギー最小化計算

ATM 変異陽性患者では PARP 阻害剤に対する有意なを実施した。

レスポンスが認められないこと、事実上全ての患者が結果および考察

治療に対する抵抗性を示すことも報告されており、臨床的に有望な PARP 阻害剤の用途を新たな視点から開いかなるタンパク質が「BRCAness」創薬標的になっ拓することが求められている。そこで、未だ詳細が明らても対応できるように、ヒトで発現している全てのタンパ

かではない「BRCAness」に関わる分子を標的とした、

ク質の WT 配列から、Alphafold2 を用いてタンパク質

革新な新規誘導阻害剤の研究開発を行うことを目的と構造を予測した。その後、申請者の小関が作成した自した。

動解析プログラム（Fortran, Bash を使用）により、各タンパク質におけるそれぞれのアミノ酸の荷電状態（プロ55（様式 20）成果報告書
トン化状態、脱プロトン化状態）及び、ジスルフィド結合
形成判定、そしてタンパク質内における水素結合ネット
ワークの最適化を実施し、生体内における最適なタン
パク初期構造を自動作成した。作成したプログラムでは、これらの初期構造を用いて分子動力学シミュレー
ション用ソフト Amber 及び Gromacs で使用できるエ
ネルギー最小化用のインプット作成までを自動で取り扱える。この作成プログラムは、後日、本申請研究とは別研究で開発中の構造解析プログラムとともに、
Gighub へのアップロードする予定である。上述した分
子動力学シミュレーション用のインプットを利用し、全タンパク質に対するエネルギー最小化計算を実施した。
申請者は、現在、ヒト全タンパク質の立体構造解析に取り掛かっている。申請者は来年度の追加利用申請をしているが、来年度では分子動力学シミュレーションを実施することで、各タンパク質の熱構造サンプリングを実施し、酵素活性部位の有無および同定と、薬物標的箇所の同定（酵素である場合は酵素部位、そうでない
場合は機能性タンパク質‐タンパク質相互作用界面の同定）をすることで、目的である BRCAness にかかわる
有効分子標的における創薬分子（低分子、中分子）の設計に取り掛かる予定である。
まとめ、今後の課題
令和３年度では、ヒトタンパク質全ての構造について
立体構造を予測し、Amber, Gromacs での計算インプットを作成してエネルギー最小化及び平衡化計算を実施した。令和４年度はこれらの結果を用いて、サンプリ
ング構造分布から、創薬ターゲットとなりえる箇所を特定し、阻害分子設計につなげる。

56

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 粗さを用いた熱・運動量輸送の非相似制御に関する研究
英文：Controlling the dissimilarity between the momentum and heat transfer by the wall roughness
利用課題責任者
桑田 祐丞
Yusuke Kuwata
大阪府立大学 工学研究科 機械系専攻
Department of Mechanical engineering
Osaka Prefecture University
URL：http://www2.me.osakafu-u.ac.jp/htlab/
邦文抄録（300 字程度）
壁面に設置された粗さが熱・運動量輸送の非相似性に与える影響を理解するために，粗面乱流熱伝達の直接数
値解析を行なった．本研究では，粗さの波長が非相似性に与える影響を調査する為に，３次元正弦波粗さを対象として，波長を変えた粗面を用いた数値解析を行った．数値解析手法に格子ボルツマン法を用い，乱流モデルを使用しない直接数値解析を行った．粗面の波長は速度場のスケーリングには影響を与えないが，温度のスケーリ
ングに大きな影響を与えることが分かった．波長の短い粗面では，スタントン数と摩擦係数との比で定義されるレ
イノルズアナロジファクタの減少が小さく，少ない流動抵抗で高い伝熱促進効果が得られることが示唆された．
英文抄録（100 words 程度）
We performed direct numerical simulation of rough wall turbulence to understand the effects of wallroughness on an analogy between the heat and momentum transfer. This study considered threedimensional sinusoidal rough surfaces with various wavelength values. We utilized the lattice
Boltzmann method for the direct numerical simulations. It is found that the wavelength does not affectthe scaling for the mean velocity but has considerable effects on the scaling for the mean temperature.
The Reynolds analogy factor, which is defined as the ratio of the doubled Stanton number to the skinfriction coefficient, is larger for the surface with a small wavelength value, suggesting that surfaces
with small wavelength roughness yield higher heat transfer rate with lower frictional resistance.

Keywords: Rough surface, Turbulence, Heat transfer, Lattice Boltzmann, Reynolds analogy
背景と目的

クタのスケーリングに関する調査を行う．

平滑面の乱流熱伝達に関しては，運動量・熱輸送の間概要に相似性が成立することが知られており，プラントル数

本研究では，少ない流動抵抗増で高い伝熱促進効が１の条件下では，スタントン数と摩擦係数との比はレ

果をもたらす理想的な粗さ構造を探求するために，イノルズ数に拠らずほぼ一定値を示す．しかし，表面粗粗面の乱流熱伝達に関する直接数値解析を実施し，さが出現すると，熱輸送に対する運動量輸送の増大が相似性を表すパラメータであるレイノルズアナロジフ

著しく，熱・運動量輸送の相似性が悪い意味で崩壊す

ァクタのスケーリングに関する調査を行う．
本研究で対象とする解析系は, 粗面を下面にもつることが知られている．つまり，粗さによる流動抵抗の増大効果は，熱伝達の増大効果を上回る．粗さを用いオープンチャネル流である．粗面の平均高さ位置た伝熱制御を考えた際，流動抵抗の増加を最小限にし𝑦 = 𝑘 から，上面 𝑦 = 𝛿 までの高さで定義される

つつ，熱伝達率を最大化することが求められるが，粗さ

有効チャネル高さ δe = 𝛿 − 𝑘 に対して，計算領域が非相似性の大きさに及ぼす影響は十分に理解されは主流方向に6δe ，スパン方向に3δe とした．入口・ていない．本研究では，少ない流動抵抗で高い伝熱促

出口境界の圧力差を一定とすることで流体を駆動

進効果をもたらす理想的な粗さ構造を探求するために，し，有効チャネル高さを代表長さとする摩擦レイノル

粗面の乱流熱伝達に関する直接数値解析を実施し，

ズ数 𝑅𝑒𝜏 = 180, 300, 600 の３つの条件で解析を行

相似性を表すパラメータであるレイノルズアナロジファった. 境界条件は, 主流方向およびスパン方向に周57

（様式 20）成果報告書
期境界条件を課し，上面は滑り境界条件を課した．
熱的境界条件として，粗面は等温条件，上面は断熱条件とし，流体には一様発熱を与えた．プラントル数
は空気を仮定し 𝑃𝑟 = 0.7 とし，浮力の影響を無視したパッシブスカラーとして温度を取り扱った．解析に用いた三次元正弦波粗面の粗面高さℎ(𝑥, 𝑧)は以下のように与えられる．
2𝜋2𝜋ℎ(𝑥, 𝑧) = 𝑘 (1 + cos ( 𝑥) cos ( 𝑧)),𝜆𝜆

図 1 波長 𝜆⁄𝑘 = 16，摩擦レイノルズ数 𝑅𝑒𝜏 =ここで, 𝜆は波長, 𝑘は振幅を表し，粗面領域は 0 <600 のケースにおける，乱流渦と温度変動の可𝑦 < 2𝑘 である．粗面の振幅 𝑘 = 𝛿 ⁄11.5 とし，波視化結果．

長𝜆と振幅𝑘との比を 𝜆⁄𝑘 = 2, 4, 8, 16 とした 4 種類の粗面で解析を行った.
熱流動の解析は 2 種類の分布関数を用いる格子
ボルツマン法によって行った．速度場は３次元２７方向多緩和時間格子ボルツマン法（1）を用いて解析し，温度場には３次元１９方向速度正規化格子ボルツマ
ン法（2）を用いる．格子ボルツマン方程式は，全て等間隔格子で離散化されており，計算格子の解像度は壁単位で 2.0 程度となるように設定した．計算格子
数は, 𝑅𝑒𝜏 = 180 のケースで 640(𝑥) × 116(y) ×320(𝑧), 𝑅𝑒𝜏 = 300 のケースで 896(𝑥) × 162(𝑦) ×
448(𝑧), 𝑅𝑒𝜏 = 600 のケースで1536(𝑥) ×277(y) × 768(𝑧)となった．プログラムコードは
CUDA Fortran によって記述されており，複数のGPU を用いた領域分割法によって並列化した．GPU間の通信は，CPU を通じて MPI 通信によって行い，分割した領域には袖領域を設けており，領域間の通信時間を隠蔽するように工夫した．
結果および考察

図１に波長 𝜆⁄𝑘 = 16 ，摩擦レイノルズ数 𝑅𝑒𝜏 =600 のケースにおける，乱流渦と温度変動を可視化した図を示す．乱流渦は速度勾配テンソルの第２不変量の等値面で可視化をしており，コンター図は温度変動を示す．乱流渦は粗面の近傍で見られ，小図 2 上図は摩擦速度で無次元化された主流平均速スケールの乱流変動による激しい温度変動が見ら

度分布，下図は摩擦温度で無次元化された平均温れる．これらの熱流動場を統計的に議論するために，度分布．壁からの距離には，Jackson モデル（３）で得

図２に主流平均速度・温度分布を示す．

られた仮想原点からの距離を用いた．

58

（様式 20）成果報告書
平均速度は，摩擦速度を用いて無次元化しており，

まとめ、今後の課題

平均温度は摩擦温度を用いて無次元化した．図２よ少ない流動抵抗増で高い伝熱促進効果をもたらす

り，粗面の平均温度・速度分布は滑面の結果と比較

理想的なパッシブ伝熱制御を確立するために，粗さして，下方に位置していることが確認できる．これは，の特徴量が運動量・熱輸送に与える影響を直接数

粗面によって壁面摩擦・熱伝達率が増加した結果で値解析によって調査した．対象とした粗面は 3 次元ある．対数領域における平均速度・温度の下方シフ

正弦波粗さとし，波長を変えたケースを複数の摩擦ト量である速度粗さ関数𝛥𝑈 + ，温度粗さ関数𝛥𝛩+ はレイノルズ数条件で解析した．直接数値解析は密度+波長が短くなるほど大きくなる．温度粗さ関数𝛥𝛩 は，

分布関数・エネルギ密度分布関数を用いた格子ボ

𝛥𝑈 + と比較して小さいが，波長𝜆⁄𝑘の影響を大きく受

ルツマン法を用いて，複数 GPU を用いた並列計算にけている様子が観察される．

よって実施した．粗面要素の波長の大きさは，速度

これら平均速度・温度の変化が流動抵抗・熱伝達場のスケーリングには影響を与えないが，温度場の率に与える影響を調べる為に，熱伝達率の無次元

スケーリングに大きな影響を与えることが分かった．

数であるスタントン数𝑆𝑡 ，流動抵抗の無次元数であ

また，波長の短い粗面は，少ない流動抵抗で高い伝る摩擦係数𝐶𝑓 を議論する．図３に，スタントン数と摩熱促進効果をもたらすことが示唆された．

擦係数との比で表されるレイノルズアナロジファクタ

参考文献

𝑅𝐴 = 2𝑆𝑡 /𝐶𝑓 と摩擦係数𝐶𝑓 との関係を示す．なお，［１］K Suga, Y Kuwata, K Takashima, R Chikasue,

図中には，レイノルズアナロジファクタ𝑅𝐴，摩擦係数

Computers & Mathematics with Applications 69 (6)𝐶𝑓 を滑面における値𝑅𝐴0 ，𝐶𝑓0 で無次元化した値を示(2015), 518-529.

す．図より，𝑅𝐴/𝑅𝐴0 は１を下回っており，粗面による［ ２ ］ K Suga, R Chikasue, Y Kuwata, International

摩擦係数の増大効果は，熱伝達率の増大効果を上journal of heat and fluid flow, 68 (2017), 225-236.

回っていることが分かった．粗面乱流熱伝達におけ［３］Jackson, P. S., Journal of fluid mechanics, 111

る速度場と温度場の非相似的なふるまいは，摩擦係(1981) 15-25.

数𝐶𝑓 が大きくなるほど顕著であり，𝑅𝐴は𝐶𝑓 の増大に従って減少することが分かる．興味深いことに，波長が短くなるほど，𝑅𝐴の減少が緩やかであることが確認される．つまり，波長の短い粗面では，より少ない
流動抵抗増で高い伝熱促進効果をもたらすと言える．

図３ レイノルズアナロジファクタ𝑅𝐴/𝑅𝐴0 と摩擦係

数𝐶𝑓 /𝐶𝑓0 との相関関係．𝑅𝐴0 ，𝐶𝑓0 は同じ摩擦レイ
ノルズ数条件における滑面の値を示す．

59

60

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

二酸化炭素の分離回収に適した分離材料の探索

Prediction of separation materials suitable for CO2 capture
藤井 達也
Tatsuya Fujii
産業技術総合研究所化学プロセス研究部門
Research Institute for Chemical Process Technology, National Institute of Advanced Industrial Scienceand Technology (AIST)
https://unit.aist.go.jp/cpt/
CO2 分離回収技術を省エネルギー化できる CO2 化学吸収液の性能を予測するモデルの構築に取り組んだ。量子化学計算により、化学吸収液に用いられる有機化合物の物性を計算し、化学吸収液に用いられる有機化合物の特徴量を獲得した。既往の化学吸収液に関する研究データを収集して、プロセスパラメータ（温度、圧力、濃度等）と有機化合物の関係をデータセット化し、既存のソフトウェアで算出した特徴量と、量子化学計算で獲得した特徴量を基に、予測モデルを構築した。モデルの性能を文献データにより検証したところ、既存の有機化合物で構成される化学吸収液については、高精度（RMSE = 0.1 未満）で CO2 吸収量を予測できた。一方で、モデルに含まれない有機化合物
を含む吸収液の性能を予測することで汎化性能を調べたところ、予測精度は低下（RMSE = 0.2 程度）した。汎化性
能の面においては予測モデルの改良やさらなる特徴量の探索が必要と考える。
Predicting method of the performance of CO2 absorbent by machine learning was investigated.
Some features for designing chemical absorption were obtained by quantum chemical calculation.
References on the performance of absorbents were collected and the dataset describing the relation
between amines, process parameters and performances. Predicting models were constructed using
features calculated by the combination of an existing descriptor calculator quantum chemical
calculation. The predicting performance was high (RMSE < 0.1) for absorbents consisting of existing
amines. On the other hand, the performance was lower (RMSE = ca. 0.2) for amine solutions that were
not included in the learning process. To improve the predicting performance, improvement of themodeling method and finding of new features should be needed.
Keywords: CO2 分離、化学吸収液、機械学習、量子化学計算、特徴量

背景と目的

吸収液設計法への発展が期待される。そこで、本研究

地球温暖化問題の解決へ向けてカーボンニュートラ

では、量子化学計算(Gaussian)と機械学習とを組み合ルを実現するためには、CO2 を省エネルギーに分離回

わせて、前記因子を入力として CO2 吸収量を予測する収する技術が必須である。分離回収プロセスにおいて、モデルを構築することを目的とした。

CO2 分離材料はプロセス全体の省エネルギー性能を左

本研究では、既往の化学吸収液に関する研究デー

右する重要な要素である。本利用課題では、CO2 分離

タを収集して、プロセスパラメータ（温度、圧力、濃度等）

回収技術を省エネルギー化できる CO2 化学吸収液のと有機化合物の関係をデータセット化した。また、プロセ

性能を、量子化学計算と機械学習の融合により予測すスパラメータ、既存の特徴量計算ソフトウェアによって得るモデルの構築に取り組んだ。化学吸収液の CO2 吸収られた特徴量に加え、量子化学計算で得られた特徴量

量は、吸収液を構成する塩基及び溶媒、塩基の組成、を組み合わせて、機械学習によって CO2 吸収量を予測

温度、CO2 分圧など数多くの因子に依存する。化学吸するモデルを構築した。構築したモデルについて、既存収液に関する従来型研究では、数多くの吸収液合成との有機化合物で構成される化学吸収液、ならびに、モ

性能試験のサイクルを回しており、開発に膨大な時間をデルに含まれない有機化合物を含む吸収液の性能の要している。CO2 吸収量の予測モデルが確立されれば、予測精度について、検証を行った。

試験数の低減など、化学吸収液の開発の加速化が期待できる。また、従来にない特徴量の発見による新規な61

（様式 20）成果報告書
概要

種類の場合と同様であることが確認された。

CO2 化学吸収液として、アミン水溶液に着目し、アミン水溶液系の既往の文献を整理して、アミンを一種類、ないし二種類含む化学吸収液について、アミンの種類、その組成、温度、圧力に対する CO2 吸収量を整理しデータセット化した。44 種類のアミンについて、4099 点のデータからなるデータセットを構築した。この中から、物理吸収の影響が大きいと考えられる高圧領域を除き、1000 kPa 未満のデータを用いてモデルを構築した(3533 点)。
図 1 アミン吸収量のパリティプロット（RDkit 特徴量を

機械学習を行うにあたり、アミンを特徴量に変換する必要がある。本研究では、既存のソフトウェアである

用いた場合）

RDkit ならびに HSPiP を用いて、特徴量を算出した。

続いて、新規アミンへの対応可能性を検討するために、

また、44 種類のアミンについて、量子化学計算を行い、

1 つのアミンをモデルの学習過程で外し、そこで学習しアミンの N 原子の電荷や、プロトン親和力等についてもたモデルでそのアミンの吸収量を予測する、ということを検討を行った。量子化学計算は、Gaussian 16 を用い、繰り返す、leave-one-group-out 法により、モデルの予

wb97xd/6-311+G(d,p)//IEF-PCM により水の分極体

測精度を検証した。その結果、44 のアミンについて、検モデル中で計算を行った。

討 し そ の RMSE の 平 均 値 を 算 出 し た と こ ろ

回帰計算は Scikit-learn のライブラリを用い、ランダ

RMSE=0.179 となった。この値は、当該アミンを含んでムフォレスト回帰により行った。予測性能の評価を行うた

学習したときの RMSE(0.079)と比べて大きな値となっめに、トレーニングデータとテストデータを 80:20 の割合ており、予測精度の低下を示している。アミン種ごとに見にランダムで分けて、トレーニングデータで構築したモると、特に環状やジアミン等の特徴を有するアミンにつデルによってテストデータの予測性を検証した。指標にいて予測精度の低下がみられており、新規特徴量やモは二乗平均平方根誤差（RMSE）を用いた。

デルの改善が予測精度の向上に寄与する可能性がある。本研究では、量子化学計算で得られた N 原子にお結果および考察

ける電荷や、プロトン親和力を特徴量として加え、予測

アミン 1 種類の系について、予測性能の検証を行っ

精度の変化を検討したところ、RMSE=0.175 とわずかた。図 1 にその結果を示す。トレーニングデータはもとよながら改善傾向も見られたため、有効な特徴量を取り入り、テストデータも良好に CO2 吸収量を再現することがれることで予測精度をさらに改善することが期待される。

分かった。RDkit 特徴量、HSPiP 特徴量を用いた場合、
RMSE=0.078, 0.077 といずれも同程度の予測性度でまとめ、今後の課題あった。ランダムフォレスト回帰においては、特徴量の重アミンの特徴量およびプロセスパラメータを利用する

要度を計算することができる。重要度を評価した結果、ことで、アミンを 1~2 種類含む水溶液系での CO2 吸収

温度、圧力等のプロセスパラメータに加えて、電荷に関量を予測可能であることを示した。一方、新規アミン種

する特徴量が CO2 吸収量の予測に重要な役割を果たが含まれると予測精度が低下することが示唆されたため、していることが明らかとなった。

モデルの汎化性能を向上することが今後の課題である。

アミン 2 種類を含む系についても、同様に予測精度

そのための方策として、モデル自体の改良に加え、新の検証を行った。その結果、アミン 1 種類の場合と同様

規特徴量を取り入れることがアプローチ法の一つとしてに、RMSE=0.079 と、高精度で CO2 吸収量を予測する考えられる。

ことができた。重要度の高い特徴量についても、アミン 162（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

光熱変換機構に関する第一原理計算
First-principles calculations on photothermal conversion mechanisms
江目 宏樹
Hiroki Gonome
山形大学 大学院 理工学研究科 機械システム工学専攻
Graduate School of Science and Engineering, Yamagata University
http://gonome-lab.yz.yamagata-u.ac.jp
太陽光は最も有力な再生エネルギーであり，これを効率良く回収・エネルギー変換を行う手段として，太陽電池や光触媒等の利用し光と物質の相互作用の理解を深める必要がある．本プロジェクトでは第一原理計算を用いて励起状態の物質の過渡的エネルギー散逸について評価を行った．イオン温度の過渡的な上昇は 1 ps 以内で起こっていることが明らかになった．α-quartz では，非結合性の電子が電子励起によって反結合性軌道を 占 有 し ， こ れ が 駆 動 力 と な っ て 結 晶 の 破 壊 を 引 き 起 こ す こ と が 示 さ れ た ． 同 様 に ， Rutile-TiO2 及 び
Anatase-TiO2 でも Fermi 準位近傍の結合電子軌道が反結合性軌道の状態に遷移することが確認された．本
プロジェクトは光と物質の相互作用の原理解明に寄与するものである．
Sunlight is the most promising renewable energy source. It is necessary to use solar cells andphotocatalysts as a means of efficiently recovering and converting this energy, and to deepen our
understanding of the interaction between light and matter. In this project, first-principles
calculations were used to evaluate the transient energy dissipation of materials in the excited state.
In α-quartz, it was shown that non-bonding electrons occupy anti-bonding orbitals due to electronic
excitation, which acts as a driving force to cause crystal breakdown. Similarly, in Rutile-TiO2 andAnatase-TiO2, it was confirmed that the bonding electron orbitals near the Fermi level transition to
the antibonding orbital state.

Keywords: 光熱変換，太陽光エネルギー，第一原理計算，分子動力学計算，時間依存密度汎関数理論

背景と目的式により計算し，逐次的にこれらのステップを繰り返し

現代の化石燃料に代わるエネルギー利用の手段と解くことで分子動力学計算を行った．また，さらに厳密して，いわゆる再生可能エネルギーを利用した技術開に励起状態の電子を取り扱うため，時間依存密度汎関発が数多く進んでいる．太陽光は最も有力な再生エネ

数理論の導入も行った．

ルギーであり，これを効率良く回収・エネルギー変換を行う手段として，太陽電池や光触媒等の利用し光と物

結果および考察

質の相互作用の理解を深める必要がある．光の励起

光と物質の相互作用を評価するため，典型的な絶縁によるキャリア輸送や熱散逸といった微視的過程は古体であるα-quartz(SiO2)，光触媒及び半導体として典的な電磁場解析では取り扱えないため，量子力学的

Rutile-TiO2，Anatase-TiO2 及び CuO を選択しそれに定式化を行う必要がある．そこで本プロジェクトではぞれ分子動力学計算を実行した．その結果，すべての第一原理計算を用いて励起状態の物質の過渡的エネ

材料でイオン温度の過渡的な上昇は 1 ps 以内で起こっ

ルギー散逸について評価を行った．

ていることが明らかになった．さらに，結合状態解析を行い電子の結合性を解析した．α-quartz では，非結概要

合性の電子が電子励起によって反結合性軌道を占有

密度汎関数理論に基づいて電子状態の解析を行い，し，これが駆動力となって結晶の破壊を引き起こすこと

さらに励起状態の電子を Fermi-Dirac 分布に従う占有が示された．同様に，Rutile-TiO2 及び Anatase-TiO2率を与えた．この時の原子の時間発展を古典運動方程でも Fermi 準位近傍の結合電子軌道が反結合性軌道63（様式 20）成果報告書の状態に遷移することが確認された．一方で，CuO は
Fermi 準位から大きく離れたエネルギー準位に結合性軌道の電子が多く存在していることが確認された．また，
Fermi 準位近傍の反結合性軌道の電子が確認されたが，これらの電子励起後の状態も同程度の反結合性を示していた．これに対して分子動力学計算を実行すると，他の材料より温度の上昇は低い値を示した．これは，電子励起による結合性軌道の電子の遷移がほとんどなく，また，遷移した電子もその前後で反結合性の差が小さく，大きな駆動力が生まれなかったためであることが示唆された．

まとめ、今後の課題
本プロジェクトでは，分子動力学計算及び結合解析
により電子と格子の相互作用に対する評価を行った．
また，光と電子の相互作用を評価するために時間依存
密度汎関数理論の導入を行った．本プロジェクトは光と物質の相互作用の原理解明に寄与するものである．

64

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 ダイラタント系のシアシックニング機構の解明に向けた理論シミュレーション
英文： A theoretical study on the mechanisms of shear thickening
利用課題責任者
山田達矢
所属
高度情報科学技術研究機構・計算科学技術部
邦文抄録（300 字程度）： 我々の研究グループでは，分子シミュレーションの方法を使用して，高分子系や粒子分散溶液系のダイラタント性（あるいはシアシック二ング）に関する理論的検討を行っている．本研究では，ポリエチレングリコール(PEG)/シリカ粒子分散水溶液系の散逸粒子動力学法(DPD)シミュレーションを行った．この系は典型的なダイラタント系で，その振る舞いは古くからよく知られているが，対応する分子レベルの現象はわかっていない．適当なモデルを用いて一連のシミュレーションを行い，構造や粘弾性に関する解析を行ってみた．
英文抄録（100 words 程度）: Shear thickening/thinning behaviors of dilatant systems of polymer solutions andnano-particle-suspended solutions have been studied in our group with molecular simulations. We focused in thisstudy on a dilatant system of polyethylene glycol solution with silica particle suspension. We analyzed molecularlevel structures and viscoelastic properties of this system through a series of dissipative particle dynamics (DPD)
simulations.
Keywords（５つ程度）: dissipative particle dynamics (DPD), ダイラタンシー，ポリエチレングリコール/シリカ粒子分
散溶液，レオロジー
背景と目的

ステップ間隔を大きく取れるため，ポリマーや生体分子

遅い剪断刺激に対しては液体として振る舞う物質が，などの大きな分子系の構造を調べるのにしばしば用い速い剪断刺激に対して固体のような抵抗力を示す場合

られる方法である．本研究でもこの理由で DPD 法を使があることが知られている（ダイラタント現象）．この種用したが，ポリマーダイラタント系への応用の前例は無の物質を衝撃緩衝材や制動装置などへ応用する研究く，ここでの調査は方法評価を兼ねる．

が行われているが，剪断速度に応じて粘度の増加する（シアシックニング）物理的メカニズムが十分に理解さ結果および考察れていないため，いまのところ材料設計の適当な指針

図１のように PEG 水溶液/シリカ粒子分散系を粗視化が確立していない．我々は，分子動力学および粗視化（DPD 粒子化）して，表１のような大小の構成のシミュレ

動力学シミュレーションを研究方法として用い，シアシッ

ーションモデルを用意した．これらモデルの構成にはダクニングの分子レベルの現象解明に取り組んでいる．

イラタント現象の観測されている実験系の物質構成[1]

本研究では，ポリエチレングリコール(PEG)/シリカ粒子

を参考にしたが，実験系とは表１脚注のような相違があ分散水溶液系のダイラタント性に関する検討を行った．

るので注意されたい．粒子間相互作用についてはあまり情報がないため，暫定的に表 2 の二種のソフトポテン

概要

シャルパラメータセットを用意した: 異種粒子が特に強この系のダイラタント性は PEG とシリカのインタープレ

い親和性や排他性を持たないもの（SP1）と，P 粒子と S

イによるものと考えられているが．構造に関する具体的粒子が親和するもの（SP2）である．

データは報告されていない．本研究では，適当なモデ
A. 計算可能なシミュレーション時間に関する検討

ルを用意して散逸粒子動力学法（DPD）シミュレーションによる解析を行ってみた．DPD 法は粗視化動力学計算

はじめに DPD 法によるダイラタント現象の直接シミュ法の一種で，通常の分子動力学計算に比較してタイム

レーションの可能性について検討した：ダイラタント現象
65

（様式 20）成果報告書の時間スケールは 10-1-10-3sec の程度で（101-103/sec の解析という側面でも（A の節で述べた事柄に加えて）検剪断速度で観測される），これは分子シミュレーション

討できる剪断速度領域が制約を受けることを象徴的にの時間スケール（タイムステップが fsec や psec オーダ示している．この問題は MD 法の計算については指摘ー）に比べるとずっと大きく，DPD 法でも実時間的に十されていたが，DPD 法でも同様であることが分かった．

分なシミュレーションができない懸念があったからであ調べた剪断速度領域では shear thinning という結果になった．また viscosity の剪断速度ゼロの極限を平衡る．
計算を実行したところ，表 2 の大小のモデル系の計系の計算から Green-Kubo 法により求めたところ，5-6h0

算速度は，スパコン 4 ノード使用でそれぞれ~2M との 値 が 得 ら れ た の で ， 0-10-6t0-1 の ど こ か で shear~20M ステップ/日であった．したがって，それぞれのモ

thickening の剪断速度領域があることが示唆される．

デ ル を 使 っ て 実 行 で き る 計 算 総 量は ， 現 実 的 に は

しかしながら計算で得られた viscosity の増減は実験系~100M と~1G ステップの程度で，シミュレーション実時に比べると小さく（実験系では値が数桁変化する），ダ

間は~1Mt0 （数十 µsec）と~10Mt0 （数百 µsec）の程度ま

イラタント現象に対応する変化かどうかは疑問である．

でと試算される（１ステップは 0.01t0, 1t0 = 25 psec, DPD

ここで得られている viscosity の数値の由来を検討す

単位については図 1 脚注参照）．非平衡動力学計算でるため，2 成分や 1 成分系の計算を行ってみた．（表 3．

は定常状態到達までにある程度大きな剪断変形を与簡易な検討のために小モデルを使用した．また統計解える必要があるが，これを仮にg=1 とすると，定常状態

析に要する計算ステップ数を少なく抑えるため，あまりを 調 べ ら れ る の は ， 小 さ く て せ い ぜ い ~104/sec や

小さくない剪断速度ということで 10-4/t0 の系を調べた）．

~103/sec の剪断速度までということになる．斯くして，ダ計算モデルからシリカを取り去ると viscosity が若干上

イラタント現象の剪断速度領域（低粘度-高粘度の遷移

昇する，また PEG を取り去ると viscosity が極端に減少すの領域）を調べるのは困難であるということが分かっる（溶媒レベルの値まで）という結果が得られていて，こた．

れは系の viscosity の殆どが PEG 由来であることを意味している．また，viscosity は PEG を寸断すると極端に減

B. 相互作用モデル SP1 を用いた DPD 計算

少するため，PEG の紐状に長いという特徴と関係したも次に，モデル系の不具合や DPD 法の計算手法に絡のであることを示している．

んだ問題などについて知見を得るため，10-2/t0 から

各剪断速度での構造を図 3 に，PEG の構造パラメー

10-6/t0 （4•108/sec から 4•104/sec）の範囲の剪断速度のタを表 4 に示した．これらのデータを全体的に見ると，シミュレーションを行い，viscosity や構造に関する一連

剪断速度の遅い非平衡系(10-6/t0)と平衡系の構造にはの解析を行ってみた．ここでは表 2 の SP1 の相互作用違いが殆ど無く，ある程度剪断速度が大きくモデルを使用した．

（10-4-10-3/t0 の辺り）なると違いが明確になるようである．

大モデル系の viscosity 解析の結果を図 2 に示した

見た目の構造については（図 3），剪断速度が 10-3/t0 以（小モデル系の結果は大モデル系と同様なので割愛）．

上で PEG 鎖が剪断速度方向に並ぶという特徴（平衡系

大きい揺らぎのために viscosity の数値をあまり正確にでは見られない）を持つことがわかる． PEG の広がり具

計算できない問題のあることが明確になった．この問題

合を表す両端距離や慣性半径（表 4）は剪断速度に応はそもそも応力計算値の揺らぎが大きいことに由来しじて比較的にはっきり変化する：慣性半径のデータを見ていて，剪断速度が小さくなるに従って viscosity（=応力ると，剪断速度が大きくなるのに伴って PEG が剪断速度

÷剪断速度）の揺らぎが深刻になっているためである

方向（x 軸方向）に引き伸ばされ，これに垂直な方向（y（応力計算値の揺らぎが剪断速度によらずにほぼ一定

および z 軸方向）には広がりを失うことがわかる．PEG のになっていることにも注意されたい）．剪断速度の小さ

結合長や結合角といった局所的な構造は，基本的に剪

い 10-6/t0 については， 100M 超のステップ数を費やし断速度に応じて値を変えない（表 4）ので，上記の変形ても大きなエラーバーが残っており， viscosity の統計

は PEG の 2 面角捩れを通じた折り畳みの具合によるも66

（様式 20）成果報告書のである．表 4 のデータを総合的に見ると，系の構造変強いられていることがわかる．この間の viscosity の数値化はどこかの剪断速度で急激に起こるというタイプのもはシリカの無い場合の値（~9.9 h0，表 3 参照，シリカののではないようである．

無い系では SP1 と SP2 の区別はないことに注意）より有以上，パラメータセット SP1 を用いた計算の結果を詳意に大きい．斯くして，速度の異なる剪断流を跨ぐよう細に調べたが，次のように実験と相容れない部分があな PEG•シリカネットワーク（図 4 でいえば縦方向の隣接

る． i) viscosity がほぼ PEG 由来でシリカとの interplay

ユニットセルに跨るもの）を持てば高 viscosity 状態を与による寄与が無い，ii) 平衡系と非平衡系で構造があえることが分かった．

まり違わず（10-6/t0 ~ 4•104/sec は実験的にはかなり大一方で ~20M ステップを過ぎると，viscosity が極端にきな速度である），iii) 違いがあっても剪断速度に応じ低下するという結果になった．これは，上記の引っ張りて連続的に変化するものである（ダイラタント転移の前のために，シリカに張り付いた PEG が徐々にはがれ，終後では急激な構造変化が起こるともの考えられている）．

には上述の剪断流を跨ぐ架橋が失われてしまったため

以上の事柄はここで使用したモデルに問題があり修正である．一旦こうなると，剪断流を跨ぐ架橋が再生するの必要があることを意味している．

ことはない（剪断流方向(横軸方向)の架橋は保持されるが）．この構造では PEG がシリカに緊密に張り付いて，

C. 相互作用モデル SP2 を用いた DPD 計算

またシリカ球が剪断流方向に整列するため，PEG の縦上記の i) の点に関係したモデルの変更として，P 粒

軸方向（剪断速度勾配方向）の空間的な広がりが 1 個

子-S 粒子間の親和性を強くしたモデル SP2（表 2）を用いのシリカ球程度しかない（図 4，表 5 の rgy_y）．このためた計算を行ってみた．ここでは前節のような大規模な検に viscosity が溶媒並みに小さい値になるというわけで討はしないで，剪断変形に対する系の応答を簡易に検ある．斯くして，1-20M ステップの高 viscosity 構造は過討するということにして，1x1x1 モデルを使用し 10-4/t0

渡的な短寿命のもので， むしろ 24M step 以降の低の剪断速度の非平衡系の計算だけを行った．

viscosity 構造が定常状態である．

計算結果を図 4 と表 5 に示した（これに対応した SP1

こうして考えると，図 4 の 1-20M ステップの構造が高による計算の結果を図 5 と表 6 に参考のために掲載しviscosity であるからと言って PEG 水溶液シリカ分散系のた）．相互作用モデル SP2 を用いた計算では，SP1 の場

ダイラタント状態に単純に対応するとは言えない．ここ

合と対照的に， PEG のほとんどのセグメントがシリカ表で行ったシミュレーションが実験に忠実なものではない

面に張り付き（図 4 と 5 の initial structure 参照）, これに（剪断速度（10-4/t0 ~ 4•106/sec）が実験系でダイラタンシ

付随してシリカ球が PEG を介して会合するのが特徴でーの観測される剪断速度（101-103 /sec）よりずっと大きある． また，PEG のセグメントの一部は，このシリカ球ク

いなどの違いがある）ことには注意する必要があるが，ラスターの間を架橋する．このために，ポリマー・シリカ

上述のような架橋には剪断速度の大小に関係なく寿命のネットワークが実効的に空間的に無限に広がった格があるはずなので，高 viscosity が維持されるメカニズム

好になっている．この構造は，PEG 水溶液シリカ分散系

が何か無くてはならない．この点で現状モデルには欠の高 viscosity 状態として ref.1 において想像されている

陥があることは明らかであり，今後の検討のためにはものに近い．

より適切なものを構築する必要がある．

この系に剪断変形を与えると，上記の架橋部が引っ張られる（架橋の上部と下部が異なる速度で横方向にまとめ、今後の課題

動くため）ことにより, viscosity の大きい状態（~20M ス

PEG/シリカ粒子分散水溶液系の DPD シミュレーショ

テップまで, h ~10-20 h0）が現れる．このタイムステップンを行った．検討の結果，以下のような知見を得た．

の間，平衡系に比べて PEG の結合角が大きい，結合偏i) DPD 法で定常化できる剪断速度系は小さくてもせ

角のエネルギーが大きい（表 5）などの特徴が認められ，

いぜい~104/sec や~103/sec の剪断速度までで，ダイラタ

剪断変形が与えられたことにより PEG が伸長の変形を

ント性の現れる剪断速度領域のシミュレーションは困難

67

（様式 20）成果報告書
である．ii) DPD 法による計算でも，MD 法の計算の場合と同様で，大きい揺らぎのために viscosity の数値を正確に計算できない問題が深刻である． iii) SP1 を SP2 の相互作用モデルを用い DPD 計算を行ってみたが，それぞれ，シリカの有無が viscosity にあまり寄与しない，定常状態が溶媒並みの viscosity しか示さない，など実験と相容れない結果を与えるものであり，今後の検討のためにはより適切なものを構築する必要がある．
参考文献:[1] M. Kamibayashi, H. Ogra, and Y. Otsubo, Journal of Colloid
and Interface Science, 321, 294–301, (2008). 表 2 のモデル構
成は，シリカ 8 vol% (9 nm 径) / PEG 0.5w%(分子量 50 万)の実
験系を参考にした．
[2] R. D. Groot and K. L. Rabone, Biophysical Journal, 81,
725-736 (2001). W 粒子の拡散係数を独自に計算し，この論文の方法に倣いt0 を算出した．
[3] Z. Shen, H. Ye, M. Kroger, and Y. Li, Phys. Chem. Chem. Phys.,
19,13294 (2017). 実際には「シリカ粒子の ダイラタント特性に関する ソフトウェア OCTA を使った 粗視化シミュレーショ
ン」 報告書の様に値を調整した．

68

（様式 20）成果報告書

図 1. PEG 水溶液/シリカ粒子分散系の DPD 粒子(P, S,W )への粗視化．
この粗視化により，エネルギー, 距離，質量，時間, viscosity の DPD
単位はそれぞれ E0 = kBT, L0 = 0.445 nm, m0 =18 amu., t0 = 24.7 psec
(1step = 10-2t0), h0 =1.1 mPa∙s になる[2].

表 1. シミュレーションモデルの構成
1x1x1
2x2x2

(-P-)1140a
4
32

(W)
87,402699,216

図 3. 2x2x2 モデルの各剪断速度での計算のスナップショット: P 粒子
(シアン)，S 粒子(赤), W 粒子は見やすさのためプロット省略.剪断速度の方向は横軸(x 軸)方向で，速度勾配の方向は縦軸(y 軸)方向である．
DPD 単位については図 1 脚注参照．

(S1139)b7
56a) P 粒子 1140 個を直鎖状に結合して PEG 分子をモデル化した(先行研究[3]の結合伸縮と偏角のポテンシャルを適用)．b) S 粒子の面心立方格子から球状(直径 ~ 9 L0)に領域を切りだしシリカクラスタを作成
した(S-S 結合に調和ポテンシャル r0 = 0.78 L0, k=16.5 E0 を適用). これ
らのモデルの構成にはダイラタント現象の観測されている実験系の物質構成[1]を参考にしたが，次のような実験系との相違がある：i)PEG 分子１本の鎖長が実験系の 1/10； ii) シリカ球の径が実験系の1/2; iii) PEG の濃度が 10 倍大きい. i)と ii)は計算系のサイズを限定的にするためのもので，iii)は，PEG 濃度が実験系のままだと，性質の違いが顕著に観測できない（純溶媒系とあまり違わない）ので，その対処である. DPD 単位については図 1 脚注参照．

表 4. 各剪断速度でのモデル系の構造パラメータ*
剪断速度 0/t0 10-6/t0 10-5/t0 10-4/t0 10-3/t0 10-2/t0
0.380 0.380
0.380 0.380 0.380 0.380l(L0)[0.016] [0.016] [0.016] [0.016] [0.016] [0.016]
140.15 140.14 140.15 140.16 140.27 140.61
q(deg.)[0.009] [0.070] [0.066] [0.066] [0.010] [0.072]
34.6
34.8
37.2
59.5
123.8 154.2
de-e(L0)[0.03] [1.37] [2.07] [2.33] [5.09] [1.08]
15.3
15.3
16.0
23.3
45.7
59.7
rgy(L0)[0.05] [0.36] [0.48] [0.62] [1.49] [0.14]
7.78
8.03
9.10
19.99 44.92 60.52
rgy_x(L0)[0.17] [0.63] [0.20] [0.82] [1.62] [2.17]
8.81
8.74
8.38
7.43
5.07
3.14
rgy_y(L0)[0.28] [0.59] [0.28] [0.37] [0.10] [0.08]
8.61
8.42
8.94
7.78
6.72
4.18
rgy_z(L0)[0.12] [0.42] [0.68] [0.22] [0.28] [0.20]

表 2.ソフトポテンシャル反発パラメ-タ(所謂 aij, E0 = kBT 単位)
SP1SP2a(P-P) a(P-W) a(P-S)
25a 26.3a
25c
a
a
25
26.3
20c

a(W-W)25a
25aa(W-S)
25c
30c

a(S-S)
35b
35ba) 先行研究モデルの値[3]を適用. b) シリカクラスタが会合しないよ
う大きめの値に設定． c) SP1 では異種 DPD 粒子が特に強い親和性
や排他性を特に持たないよう設定；一方 SP1 では P 粒子と S 粒子が親和的になるよう設定．

*各エントリ: PEG 結合長と結合角（l,q）, PEG 両端距離（de-e）, PEG 慣性
半径とその x-z 成分（rgy, rgy_x など，剪断速度は x 軸方向,速度勾配は y軸方向）. [ ]内数値は標準誤差．数値は DPD 単位で表現(DPD 単位については図 1 脚注参照)．

図 4. SP2 相互作用モデルを用いた計算における各ステップでの構造．
P 粒子(シアン，紫，緑. ピンクで 4 個の PEG 分子を区別)，S 粒子(赤),W 粒子は見やすさのためプロット省略. 剪断速度（10-4/t0 で計算）の方向は横(x 軸)方向で，速度勾配の方向は縦(y 軸)方向．初期構造は平衡系のスナップショットから取った．数値は viscosity の 4M ステップ
区間平均と標準誤差（DPD 単位については図 1 脚注参照）．

図 2. 2x2x2 モデル系の各剪断速度での viscosity 解析の結果： 記載の数値は平均値,標準誤差(角括弧)と計算ステップ数（丸括弧)．h(+0)
は Green-Kubo 法に基づき平衡系の計算データから解析．DPD 単位については図 1 脚注参照．

表 3. 構成の異なるモデル系の viscosity (h0 単位)*
P+W+S P+W
W+S
W
P’+W+S
7.70
9.94
0.88
0.55
0.98
[0.46] [0.63] [0.15] [0.17] [0.18]*1x1x1 モデルを使用し剪断速度 10-4/t0 で計算. [ ]は標準誤差.
P’: ポリマ-結合を寸断し２粒子分子にした．h0 =1.1 mPa∙s.

69

（様式 20）成果報告書
表 5. SP2 相互作用モデルによる計算結果（10-4/t0 で計算）

表 6. 前段の続き

step
h(h0)
ref.*
10.1
0-4M
17.5
4M-8M 21.4
8M-12M 11.3
12M-16M 9.11
16M-20M 13.8
20M-24M 1.31
24M-28M -0.29

stepde-e(L0)
ref.*
42.1
0-4M
61.5
4M-8M
83.9
8M-12M
89.4
12M-16M 66.6
16M-20M 74.5
20M-24M 80.8
24M-28M 58.5l(L0)0.38000.38000.38000.38000.38000.38000.38000.3800q(deg.) Ebond(E0) Eangle(E0) Enb(E0)140.18 15862.5
140.93 15856.1
141.16 15855.8
141.10 15853.8
141.12 15854.2
141.15 15853.5
141.04 15855.0
141.03 15853.1

2322.29 468867
2325.71 468885
2332.54 468911
2331.95 468886
2331.94 468852
2332.72 468852
2329.27 468830
2328.92 468786

rgy_x(L0)
10.7
17.0
27.4
28.7
22.2
25.1
25.5
19.8

rgy_y(L0)
9.1
12.5
9.9
7.9
7.1
8.3
8.3
6.6

rgy_z(L0)
8.6
6.9
6.6
8.6
8.3
7.9
8.5
6.7

*ref.: 平衡系(剪断速度 0)の結果．各エントリ: viscosity(h), PEG 結合

表 5. 前段の続き
stepde-e(L0)
ref.*
28.8
0-4M
52.3
4M-8M
97.1
8M-12M
90.6
12M-16M 88.8
16M-20M 88.3
20M-24M 68.7
24M-28M 51.9

rgy(L0)
17.4
23.2
30.6
31.5
25.1
28.1
28.9
22.5

rgy(L0)
13.5
21.4
34.0
31.5
31.0
30.5
25.9
22.2

rgy_x(L0)
4.2
16.0
30.5
28.5
27.4
27.3
22.7
18.9

rgy_y(L0)
9.6
8.7
7.4
5.8
5.4
5.4
3.9
4.1

長と結合角(l,q), ポテンシャルエネルギー成分(Exxx,結合伸縮・偏角とrgy_z(L0)
7.0
7.6
6.8
6.3
6.9
6.3
6.7
6.8

非結合相互作用), PEG 両端距離(de-e), PEG 慣性半径とその成分(rgy,剪断速度は x 軸方向,速度勾配は y 軸方向). 数値は DPD 単位（図 1
脚注参照）で表現.

*ref.: 平衡系(剪断速度 0)の結果．各エントリ: viscosity(h), PEG 結合長と結合角(l,q), ポテンシャルエネルギー成分(Exxx,結合伸縮・偏角と非結合相互作用), PEG 両端距離(de-e), PEG 慣性半径とその成分(rgy,剪断速度は x 軸方向,速度勾配は y 軸方向). 数値は DPD 単位（図 1
脚注参照）で表現.

図 5. SP1 相互作用モデルを用いた計算における各ステップでの構造．
P 粒子(シアン，紫，緑. ピンクで 4 個の PEG 分子を区別)，S 粒子(赤),W 粒子は見やすさのためプロット省略. 剪断速度（10-4/t0 で計算）の方向は横(x 軸)方向で，速度勾配の方向は縦(y 軸)方向．初期構造は平衡系のスナップショットから取った．数値は viscosity の 4M ステップ
区間平均と標準誤差（DPD 単位については図 1 脚注参照）．

表 6. SP1 相互作用モデルによる計算結果（10-4/t0 で計算）
steph(h0) l(L0) q(deg.)ref.*~2 0.3800 140.12
0-4M
10.8 0.3800 140.18
4M-8M 17.6 0.3800 140.21
8M-12M 10.2 0.3800 140.16
12M-16M 6.7 0.3800 140.11
16M-20M 9.5 0.3800 140.15
20M-24M 8.5 0.3800 140.18
24M-28M 7.4 0.3800 140.17

Ebond(E0) Eangle(E0) Enb(E0)16080.1 2318.38 469798
16079.1 2318.78 469797
16080.0 2322.17 469785
16079.4 2321.52 469793
16082.0 2319.05 469802
16080.2 2320.31 469798
16079.9 2320.86 469794
16080.6 2319.96 469796

70

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 二酸化炭素の分離回収におけるゼオライト構造の最適化
英文：Optimization of zeolite structure for CO2 separation
利用課題責任者 池田歩
Ayumi Ikeda
所属 産業技術総合研究所
National Institute for Advanced Industrial Science and Technology (AIST)
https://www.aist.go.jp/

邦文抄録（300 字程度）
二酸化炭素吸着に適したゼオライト構造を探索するため、二酸化炭素の分子径に近い細孔径をもつゼオライト
について、第一原理計算に基づく密度汎関数理論 (DFT) 計算を利用して、二酸化炭素の吸着エネルギーを算出した。まず、ゼオライトの組成 (Si/Al 比) ごとに、最安定エネルギーを計算し、各組成のゼオライト中に二酸化炭素を存在させたときのエネルギーを計算した。この結果、算出したチャバザイト型ゼオライトの CO2 吸着
エネルギーは、実験値と良好に一致することを確認した。これより、DFT 計算を利用することで、分子動力学法やモンテカルロ法と異なり、ノンパラメータでゼオライトに対する CO2 吸着エネルギーが算出できることを示した。
英文抄録（100 words 程度）
Searching for zeolite structures suitable for carbon dioxide adsorption, we used DFT calculations
based on the first-principles calculations to calculate the adsorption energy of carbon dioxide forzeolites with pore diameters close to the molecular diameter of carbon dioxide. First, thestructure-optimized energy was calculated for each zeolite composition, and then the energy ofcarbon dioxide in zeolite structure was calculated. As a result, the calculated CO2 adsorption energyof chabazite zeolite was confirmed to be in good agreement with the experimental value. This showsthat the DFT calculation can calculate the CO2 adsorption energy for zeolite with non-parameters.

Keywords: DFT calculation, CO2 adsorbent, zeolite

背景と目的料の設計指針を確立でき、吸着材開発の加速化を期

カーボンニュートラルの実現には、省エネルギーな待している。

CO2 分離回収技術が必要である。分離回収プロセスで使用される CO2 分離材料は、CO2 分離回収プロセスの概要

省エネルギー性能に大きな影響を及ぼす。ここでは、

本研究では、TSUBAME を利用して第一原理計算

CO2 分離材料としてゼオライトに注目する。ゼオライトに基づく DFT 計算（RSDFT）を実行できる環境を整え

は、直径が 1 ナノメートル未満の細孔を有する結晶性た。数種のゼオライトについて、アルミニウム含有量を無機材料であり、その CO2 吸着性能は細孔容積と組成変えた際の最も安定な構造とそのエネルギーCO2 をゼによって決定される。細孔容積は骨格構造に依存し、オライト骨格内に配置したときのエネルギーを計算し、組成は骨格内のアルミニウム含有量や対カチオンによ

その差から吸着エネルギーを算出した。この結果、って変化することが知られている。 本研究では、第一0.38 nm の細孔をもつチャバザイト型ゼオライトの CO2原理計算に基づく DFT 計算（RSDFT）を利用して、ゼ吸着エネルギーが実験値と良好に一致することを示し

オライト中の Al やカチオンの安定位置の特定、CO2 吸

た。

着エネルギーの算出を行う。本研究により、DFT 計算
によって、ゼオライトの構造及び組成が CO2 との相互結果および考察

材用に及ぼす影響を解明できれば、固体 CO2 吸着材

はじめに、8 員環をもつチャバザイト型ゼオライトにつ71

（様式 20）成果報告書
いて、アルミニウム(Al)原子とナトリウムイオンの安定る影響を検討した (表 1)。このとき、CO2 は二重 6 員環

位置を決定した。図 1 は、チャバザイト型ゼオライトのに対して垂直に配置して計算した。表 1 より、二重 6 員骨格構造を示す。チャバザイト型ゼオライトのユニット

環の上面にナトリウムイオンが存在する場合（A,B）、ナ

セルの組成は NaxAlxSi36-xO72 とし、Al 原子数 x=1, 2

トリウムイオンが二重 6 員環内にある場合 (C) と比較(Si/Al 比=35, 17) で計算した。Si/Al 比＝35 (ユニット

して CO2 吸着エネルギーは高いことが分かった。ここで、セル内の Al 原子 1 個) のチャバザイト型ゼオライトに

実験値の CO2 吸着熱は 37.4 kJ mol-1 だったことから、

ついて、Al 原子を二重 6 員環と 4 員環の交点に配置し現実のチャバザイト型ゼオライト膜中のナトリウムイオ

たとき、ナトリウムイオンが 4 員環より 6 員環近傍にあンは A や B の構造に近い可能性があると推察した。

る方が-0.53 eV 安定であった。Al 原子とナトリウムイオンの距離に着目すると、ナトリウムイオンを 4 員環近傍

表 1 Al 原子と Na イオンの位置が CO2 吸着エネルギに配置した場合、距離が近いほど構造最適化エネルギ

ーに与える影響 (Si/Al 比=35, 黄: Na イオン, 赤: Al)ーは低くなり、安定であることが分かった。一方、ナトリ
ウムイオンを二重 6 員環近傍に配置した場合、Al 原子とナトリウムイオンの距離が構造最適化エネルギーに有意な差を与えなかった。以後、二重 6 員環近傍にナトリウムイオンを配置した場合について考える。
Si/Al 比=17 (ユニットセル内の Al 原子 2 個) のとき、次に、Si/Al 比の異なるチャバザイト型ゼオライトについナトリウムイオンが 6 員環近傍かつ各イオンが最も離て CO2 の吸着熱を DFT 計算により求めたところれた位置にあるとき、もっともエネルギーが低くなった。

Si/Al=3~∞において実験結果とよく一致することが明確になった。

まとめ、今後の課題
DFT 計算を利用することで、力学パラメータ無しで、チャバサイト型ゼオライトに関する CO2 吸着エネルギー
および拡散係数を算出できるようになった。今後、吸着等温線を描画するため、圧力依存性と温度依存性を表
図 1 チャバザイト型ゼオライトの構造[1]現することが課題である。

次に、チャバザイト型ゼオライトに対する CO2 の吸着

Reference[1] Ch. Baerlocher and L.B. McCusker, Database

エネルギーを算出した。式(1)に、CO2 の吸着エネルギ
ーつまり、吸着熱 ΔHads と DFT 計算から得られるエネ

of Zeolite Structures:

ルギーの関係を示す。

http://www.iza-structure.org/databases/(

H ads = Ezeolite+CO2 − Ezeolite + ECO2)

・・・(1)

ここで、Ezeolite は構造最適化したゼオライト構造のエネ
ルギー、Ezeolite+CO2 は CO2 を配置して構造最適化した
ゼオライト構造のエネルギー、ECO2 は CO2 分子がもつエネルギーを表す。
Si/Al 比=35 (ユニットセル内の Al 原子 1 個) のとき、Al 原子とナトリウムイオンの位置が CO2 吸着熱に与え72

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

利用課題名 臨床情報統合データベースの機械学習解析
英文：Machine Learning for Integrated Database of Clinical Information
利用課題責任者 浅井 聰
Satoshi Asai所属 日本大学 医学部 生体機能医学系 薬理学分野
Department of Pharmacology, School of Medicine, Nihon University
https://nu-pharmacology.com
本課題は医療データベースから生活習慣病に関する知見を抽出するためにノンパラメトリックベイズ法の効率的な大規模アルゴリムを実装することを目的とする。特にノンパラメトリックベイズ法においてスパース性と潜在変数を扱えるモンテカルロアルゴリズムを構築し、データの背後にある潜在的医学生物学的メカニズムを抽出しつつこれに最も影響を与えるリスク要因を明らかにしようとする。この目的で、申請者のグループが開発した新奇な手法を GPU アクセラレーターを利用して大規模データベースに適用可能にする。本年度は既存のアルゴリズムを GPU にオフロードし高速化することに成功し、これを改変することで提案法のコードを開発した。
We aim at extracting information about biophysiological mechanisms and key factors for the onset
and progression of metabolic diseases from a large volume of electric medical records. For this purpose,
we develop an efficient, sparse nonparametric Bayesian algorithm with latent variables based onMonte-Carlo sampling. In the fiscal year 2021, we succeeded to develop a code that drasticallyaccelerates a conventional algorithm by using GPU accelerators, and then, we developed a code for our
new method by modifying this code.

Keywords: electric medical record, nonparametric Bayes, Monte-Carlo algorithm, GPU computing,metabolic diseases
背景と目的て推定できることも必要だ。これらを変数間の非線形関

近年、電子カルテデータなどを含む大規模な医療デ係も考慮しながら行うことは，潜在変数付きノンパラメト

ータベースから薬剤の効果・疾患リスクなどの医学的リックベイズモデルでスパースな変数選択を行う問題に知見を機械学習・人工知能技術を用いて抽出する試み

帰着される(次ページ図参照)。上記のような潜在変数

が注目を集めている。実際に、世界各国の大学病院や付きのノンパラメトリックベイズの枠組みにおいて、医療

地域中核病院の電子カルテデータを用いた解析が行

データベースのような大規模データ上でスパースかつ

われはじめており[1,2]、さらに適切な倫理的な枠組み

効率的な推定に成功した研究はまだ存在しない。先行のもと、電子カルテデータに患者から採取した遺伝•生

研究のほとんどで潜在変数のモデリングは巧妙に避け

化学的測定データを組み合わせ、これを患者個別医療られている。実際、医療データに限らずノンパラメトリッに役立てようとする試みがはじまっている。

クベイズモデルにおける潜在変数の推定やスパース変

しかしながら、近年の機械学習•人工知能技術の発数選択は、機械学習一般の問題として依然として困難展にもかかわらず、上記のような大規模医療データかな問題の一つであり、効率的な解法が模索されている

ら有益な医学的知見を抽出するためには解決すべき[3]。申請者のグループは、下にも述べるように、この問

技術的課題が複数存在する。例えば、医療データベー

題を解決しうる新奇なアルゴリズムと一定の理論的保スの解析では薬剤投与・検査値推移・疾患発症など数証を得たので、本課題において TSUBAME を利用し千以上の予測変数から疾患進行を予測しようとするが、た大規模実装を行いこの手法の有用性を示そうとする。

実際に疾患進行に影響する変数は少数であり、変数選

本年度はこの実装に一定の進展があったので報告す択を効率的に行う必要がある。また臨床医が判断するる。

重症度のような潜在変数の情報をアルゴリズムによっ73

（様式 20）成果報告書可能となった。この従来法の実装は申請者らの提案
概要

法との比較において役立つとともに、この実装の改

上記「背景と目的」に述べたように、大規模医療データ

変によって提案法を実装していくための実験コードとベースから医学的知見を引き出すには潜在変数付き

して役立つと期待される。実際、この結果から申請のノンパラメトリックベイズモデルにおいて効率的な変者らは提案法の場合の大規模実装の見通しを得る

数のスパース選択を可能にする必要がある。申請者のとともに、実装されたコードを改変しながら提案法のグループではこれを可能にするモンテカルロアルゴリズ

コードの開発も進めることができた。

ムを設計し、その理論的性能保証を得た。そこで実際に大規模医療データベースに適用するための大規模

まとめ、今後の課題実装を開発し、その性能を実証しようとする。特にTSUBAME での課題に取り組み始めて 4 ヶ月であTSUBAME の複数 GPU 環境を利用することで大規模り課題解決にさらに時間を要するが、もう数ヶ月で提データでも現実的な時間内に推論が可能であることを案法の実装が完了する見通しである。完了後はベン示し, さらにアルゴリズム面でも従来法との比較におい

チマーク用データを用いてアルゴリズムの性能を従て提案手法の優位性を示す. その後, 倫理委員会の来法と比較し、提案法の有用性を示し論文発表する。

承認を経た上で実際の匿名化データにアルゴリズムを

また、倫理委員会の承認を経て実際の医療データに適用し医学的有用性を示す.

アルゴリズムを適用していく。

結果および考察

参考文献

提案法を実装する準備として，非ベイズの設定に[1]S. Lee and H.-S. Kim J Lipid, Atheroscler.

おけるマルチカーネル学習の従来法を(2021) 10(3):282-290

TSUBAME3.0 上で大規模実装し, 動作確認を得[2]M. Chowdhury et al., Front. Psychiatry

た。電子カルテデータで想定される百万サンプル×(2021) 738466

数百〜数千予測変数での計算を C++と Openacc に[3]M. Gönen, Proceedings of ICML (2012)より，全て GPU にオフロードすることで高速計算が74

（様式 20）成果報告書

TSUBAME 共同利用 令和３年度 学術利用 成果報告書

CFD 解析を用いた卓球ボールの空力特性の計算
Calculation of aerodynamic characteristics of table tennis ball using CFD analysis
伊藤建一
Kenichi Ito
新潟工科大学Niigata Institute of Technology
https://www.niit.ac.jp/
申請者は，有限要素法ソフト COMSOL を用いて，卓球ボール周りの流れを CFD 解析することにより，高精度に空力特性計算したいと考えている．本申請では，TSUBAME の実行環境整備及びアプリケーションの利用可能性について調査を行なった．その結果，COMSOL の遠隔実行を成功させ，さらに実行速度がどの程度改善されるか確認した．
The applicant would like to calculate the aerodynamic characteristics of a table tennis ball with highaccuracy by CFD analysis of the flow around the ball using the finite element method software
COMSOL. In this application, we investigated the execution environment of TSUBAME and thepossibility of using the application. As a result, we succeeded in remotely executing COMSOL, and also
confirmed how much the execution speed could be improved.

Keywords: table tennis, aerodynamic, CFD, COMSOL

2.

背景と目的
卓球ボールは小型軽量であり，ラケットに貼付されて研究室で作成した COMSOL の流体サンプ

ルモデルを用いて複数のノードで実行し，実行速いるゴム製のラバーによってボールに回転をかけやす

度がどの程度改善されるか確認した．

いため，ボールに働く空気力の影響によって飛翔する
軌道が大きく変化するという特徴がある．卓球の競技

概要
TSUBAME の利用開始日は令和 4 年 1 月 28 日で力向上には，空気力の影響によってどのように軌道があり，利用可能期間は 2 か月程度であった．また，申請変化するかを把握しておくことが求められる．
申請者は，CFD 解析によって求められた空力特性

者は TSUBAME のようなスパコンは未経験であったに物理運動特性を考慮した飛翔軌道シミュレーターのので，今回の申請の主目的は，TSUBAME 上での開発を試み ている ．本研究では，有限要素法ソフト

COMSOL の実行環境整備であった．なお，COMSOLCOMSOL を用いて，卓球ボール周りの流れを CFD 解ライセンスはトライアルのため有効期間が 1 か月であっ析することにより，高精度に空力特性計算する予定でた．この期間内で，申請者は上記目的を達成することある．この CFD 解析にはマシンパワーが必要な計算をができた．

多くのパターンで実行する必要がある．本申請では，また，計算ノード数を増やすことによる計算時間の改

TSUBAME の実行環境整備及びアプリケーションの善についても 2 種類のモデル（モデル A とモデル B）を利用可能性について調査を行ない，以下の成果を得用いて簡易的に確認した．モデル A は，研究室のパソ

た．

コン（CPU：Core i7-7820X 3.6GHz）で 802 秒と短時
1.

研究室のパソコンで COMSOL ライセンス間で計算が終了するもので，TSUBAME のノード数をサーバーを起動し，TSUBAME 上にインストール

1, 4, 6, 8, 16 に変更して計算時間を確認した．モデル

された COMSOL を遠隔から実行することに成功

B は，研究室のパソコンで 1376 分と比較的長時間で計した．

算が終了するもので，TSUBAME のノード数を 1, 4, 8,

75

（様式 20）成果報告書
10, 12, 16 に変更して計算時間を確認した（ただし，ノード数 16 ではエラーにより計算が終了しなかった．）．

結果および考察
図 1 と図 2 に各モデルの計算時間を示す．どちらのモデルにおいても，計算ノード数 8 個の時が最も計算時間が短くなった．計算時間は，計算ノード 1 個と比較して，それぞれ約 42%と約 29%まで減少した．また，研究室のパソコンと比較すると，それぞれ約 31%と約15%まで減少した．なお，計算結果は，計算ノードを増やしてもほぼ同じであった．
700

計算時間 [s]

600

500
400

300
200
100
0
0

5

10

15

20

計算ノード数

図 1 モデル A の計算時間
800
計算時間 [min]

700
600
500
400
300
200
100
0
0

5

10

15

20

計算ノード数

図 2 モデル B の計算時間
まとめ、今後の課題
今年度の利用期間は 2 か月と短かったため，当初から COMSOL の実行環境整備及びアプリケーションの利用可能性の調査を目的とした．来年度以降は，今年度の成果にもとづいて，実際に COMSOL を用いて空力特性を計算したいと考えている．

76

東京工業大学 TSUBAME 共同利用
令和３年度利用終了課題 利用成果報告書集

発行 ： 令和５年 ５月
国立大学法人 東京工業大学 学術国際情報センター 共同利用推進室
住所 ： 〒152-8550 東京都目黒区大岡山 2-12-1 E2-6
E-mail : kyoyo@gsic.titech.ac.jp
URL : https://www.gsic.titech.ac.jp/tsubame
本書に記載の記事・写真等の二次利用を禁じます。これらの情報は著作権法上認められた「私的利用」または「引用」の条件をみたした場合を除いて、著作権者に無断で転載、複製、
放送、公衆送信、翻訳、販売、賃与等の利用を禁じます。