学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

jh200004-NAH

分子性結晶における水素ダイナミクスと同位体効果の起源解明
立川 仁典（横浜市立大学）
概要
本研究では、経路積分分子動力学(PIMD)プログラムと、固体や結晶の電子状態計算を可能とする第一原理計算法を効率的にリンクさせるためのプログラム開発を行った。このような HPCI に最適な並列プログラムを実装・使用することにより、水素を重水素に置換することで生じる大きなH/D 同位体効果の発現機構の解明を目指した。特に、本研究では有機性結晶である四角酸結晶を対象とした。四角酸結晶(H2SQ)は二次元の水素結合ネットワークを形成しており、現在、注目を集めている。H2SQ の水素結合構造に着目し、H2SQ と D2SQ における同位体効果について解析した。シミュレーションによる解析の結果、水素結合長は D2SQ の方がより大きく変化する傾向にあることを見出した。

1. 共同研究に関する情報
2. 研究の目的と意義(1) 共同研究を実施した拠点名
九州大学・情報基盤研究開発センター

研究の目的
経路積分分子動力学(PIMD)プログラムと、(2) 共同研究分野

固体や結晶の電子状態計算を可能とする第超大規模数値計算系応用分野

一原理計算手法を効率的にリンクさせることにより、HPCI に最適な解析的並列プログラムを実装すること 1 つ目の目的とした。次(3) 参加研究者の役割分担
全体統括：に開発したプログラムを用いることにより、

立川仁典 (横浜市立大学)

分子性結晶中の水素結合構造が引き起こす

PIMD プログラム開発とチューニング：

構造相転移や重水素に置換することで生じ小林理 (横浜市立大学)

る大きな H/D 同位体効果の発現機構を明ら

高木牧人 (横浜市立大学)

かにすることを 2 つ目の目的として研究を行兼松佑典 (広島大学)った。

PIMD 計算実行と解析：
石元孝佳 (広島大学)

研究の意義

David Rivera Samuel Rocabaro『水素結合』は、自然界で最も単純な結合(横浜市立大学)様式であるにもかかわらず、生体物質や水なAlbert Iskandarov

ど我々の存在と深く関わっている。ところが(横浜市立大学)

この水素結合の結合原理は未だに不明な部

坂上弘輝 (横浜市立大学)

分が多く、これまでにも不可思議な現象が数多く発見されている。特に分子性結晶中の水1

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

素・重水素置換による構造相転移温度の著しい変化や水素結合構造の変化における H/D

3.

同位体効果の解釈には、プロトンの量子力学

当拠点公募型研究として実施した意義
本研究では、我々がこれまで開発してきた的取り扱いが重要となる。

独自理論手法をさらに発展させ、高精度化・

例えば四角酸(H2SQ: C4H2O4)結晶(図 1)で大規模化計算を可能とする on-the-fly PIMD

は常温時には強誘電体であるが、温度上昇にを開発することを目的としている。

より反強誘電体へと相転移するこことが知PIMD では、従来の結晶系に対する第一原られている。H2SQ の二次元水素結合ネット

理計算では直接取り込むことの困難な水素

ワーク中の水素原子は、片方の酸素原子に結の熱・量子揺らぎの効果を取り扱うことを可合した構造から相転移に伴って二つの酸素能とする手法である(図 3)。また、この手法を原子の中央に位置した構造に変化すること用いることにより、重水化された場合の量子が実験から見出されている(図 2)。また、結晶効果がどのように変化するかについても、定を構成する水素を重水素に置換した D2SQ で量的に解析が可能となる。

は、相転移温度が約 150K も上昇することがPIMD では、図 4 に示すように、原子核の知られている。しかしながら、実験のみでは量子性を古典粒子の集まり(ビーズ)として実水素原子の振る舞いを直接解析することが現することで、各ビーズの第一原理計算を独困難である。そのため、この相転移温度に対立に取り扱うことが可能となり、高い並列化する H/D 同位体効果の機構解明に至ってお効率が実現できる。このように PIMD は並列らず、理論計算を用いた解析が望まれている。

計算機システムの性能を最大限に引き出すことが可能である。本研究で開発を行った階層的並列プログラム on-the-fly PIMD では、PIMD 部分だけでなく、第一原理計算部分にも並列処理を採用することで高い並列化効率を維持し、大規模計算を実現することが可能である。

図 1 四角酸結晶 (C4H2O4)n (H-SQ)

図 3．第一原理経路積分法の概略（複数の古典粒子の集まりで、量子揺らぎを表現）

図 2 H-SQ は相転移で水素原子の位置が変わる。重水素化すると、相転移温度が 150℃も上昇する。このメカニズムは現在のところ不明。
2

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

ステムにおいて最大 32 ノード並列計算を実行した経験を有していた。加えて、PIMD とVASP を連結させたプログラムを ITO サブシ
ステムの数ノードという小規模サイズではあるが、プログラム実装に着手していた。そのため、開発プログラムを 32 ノードで実行する際には、大きな問題は生じることなく研究計画どおりにプログラムの開発および簡単なテストを実施することが出来た。プログラムの開図 4．第一原理経路積分法のアルゴリズム

発後は PIMD-VASP プログラムにおいて、既存（経路積分と第一原理計算によるデータでの結果検証と共に、並列化効率など階層的並列構造）の向上に向けたプログラムの改良に取り組んだ。
中期(令和 2 年 8～9 月)には、大規模計算のこれまでに我々は、九州大学計算機システ実行に向けたテスト計算を中心に取り組んだ。

ム ITO A において、第一原理計算プログラム

まずは計算モデルの単純な H2SQ を例に大規VASP (Vienna ab-initio simulation package)を用模計算の安定的な動作確認と適切なシミュレいた大規模第一原理計算を実行した実績をーション条件の検証を行った。このテスト計有している。既に VASP は九州大学計算機シ算期間中に計算結果の不備等があればプログステムに導入されているが、本研究では、さラムの再確認を行う予定であったが、大きならなる高効率・大規模計算を実現するために、コンパイル時の最適化オプションの選別等、問題が生じることなく計画通りにシミュレーション・プログラムを用いた検証を実施する

ITO A 計 算 機 環 境 に 最も 適 し た 条 件 で のことが出来た。

VASP の実行環境を整えている。よって、速後期(令和 2 年 10 月～令和 3 年 3 月)では、

やかに on-the-fly PIMD のテスト計算および具体的な応用計算に取り組んだ。大規模シミ応用計算に取り組むためにも、九州大学計算ュレーションによって、H2SQ/D2SQ では大き機環境を活用することが不可欠であった。以な構造相転移温度の変化と水素/重水素の量上のような観点から、本研究では、九州大学

子効果・温度効果の関係を調べるために、様々情報基盤研究開発センターにおいて、公募型な条件下でシミュレーション計算を実施した。

共同研究を実施した。

シミュレーション研究を実施するために(最長 1 週間の計算時間制限のため)、1 週間ご

4. 今年度の研究成果の詳細とにリスタートして計算機資源を有効活用し計算機利用上の工夫た。また、並列化に関する検討も事前に行い、研究機関の前期(令和 4 月～7 月)には、開発十分な並列性能が得られるようチューニングプログラムの実装および大規模計算に対するした。しかしながら、VASP の計算時に、各ノテスト計算を実施した。我々は、九州大学情報ードで初期構造に違いから収束時間にばらつ

基 盤研 究開発 セン ターの ITO A にお いてきが見られることもあり、VASP の収束性の向PIMD プログラムを移植した経験を有してい上のための計算手法の適切な選択などが必要る。また、結晶系の電子状態計算に使用予定でであった。本研究では、このような工夫を行うある VASP プログラムに関してもすでに同シことによってシミュレーションを実施した。
3

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

これまでのシミュレーションではクラスタ算環境が必要であるため、具体的に計算を進ーモデルのみの解析が多く、結晶全体の物性めるに至っていなかった。そこで、本研究でと水素結合の関係性が明らかではなかった。

は、九州大学・情報基盤研究開発センターの本研究により、結晶構造中の水素結合様式を大型スパーコンピュータを用いてシミュレー直接取り扱うことが可能とするための基礎がションの立場より現象解明を目指し、そのた構築できたと考えている。特に、本研究を通しめに必要となるプログラムの整備を行った。て、特異な伝導率の出現機構と水素ダイナミ

次に本研究で見出したサイエンスの成果につクスの関係性を解明するための道筋についていて記載する。

多くの事柄が明らかになってきた。

本研究では、分子性結晶である四角酸(H2SQ)を計算対象とし、結晶中に存在する低障壁水素結合により発現する特異な構造相転研究成果

移温度における水素ダイナミクス、および顕これまでに経路積分分子動力学(PIMD)プ著な H/D 同位体効果の機構解明に取り組んだ。

ログラムの開発および実装に取り組み、電子四角酸結晶は図 5 に示すような有機性結晶と核の全自由度を量子力学的に取り扱ったであり、二次元の水素結合ネットワークを形高精度計算を実施している。この手法は原子成している。本研究では、特に H2SQ の水素結核の量子揺らぎや温度効果を直接考慮する合構造に着目し、H2SQ と D2SQ における同位ことのできる有効な手法である。特に、水素体効果について解析を目指した。特に、本研の量子効果は低障壁型水素結合(LBHB)で極究では、酸素・酸素間距離、酸素・水素間距離めて重要であることが指摘されている。例えに注目して研究を実施した。以下、結果につば、森ら(東大)は分子性結晶である H3(Cat-いて記載していく。

EDT-TTF)2 において、世界最高温度での伝導度や重水素置換体での特異な絶縁体相転移を見出している(Nat. Comm., 4, 1344 (2013))
が、分子論的な発現機構の解明には至っていない。この問題に対して、これまでに PIMDプログラムを用いて、H3(Cat-EDT-TTF)2 のク

RO1Hc

ラスターモデルに対して開発手法を適用し、
b量子化学計算を行うことに成功している。こaRO2H

ROO
δH(D) = RO1H - RO2Hのような計算結果から水素の量子効果や-型分散力の寄与の重要性を見出したが、伝導

図 5．四角酸結晶と水素結合部分パラメータ性や相転移についてはより大規模な計算が必要となることも同時に判明した。

多数のシミュレーションを行った結果、水これまでの PIMD による量子化学計算には、素と重水素の違いを考察していくためには、孤立分子やクラスターなどの取り扱いを得意水素（もしくは重水素）のみに注目するだけとする Gaussian などのプログラムを用いてきでは十分ではなく、周辺に存在する酸素等のた。しかし、分子性結晶を取り扱うには VASP原子にも着目していく必要があることが明などの第一原理計算プログラムとの連結が必らかになった。特に、水素と酸素の振る舞い要である。しかしながら、数百原子を含む分について温度を変化させたときに、シミュレ子性結晶モデルの第一原理計算には大規模計ーション中でどのように変化していくかに4

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

ついて着目していくことが重要であることが判明した。

a) H-SQ

図 6 には、300K における酸素間の分布をまとめた。ここで、図 6a には H-SQ を、図 6b には重水素化された D-SQ の計算結果を示す。
a)H-SQb) D-SQb) D-SQ

図 7 600K のときの酸素間距離の 1 次元分布。
図中の a) は H-SQのシミュレーションから得られた結果であり、b) は重水素化された D-SQ のシミュレーション結果から得られたものである。
図 6 300K のときの酸素間距離の 1 次元分布。
図中の a) は H-SQ

図 6、図 7 から分かるように、通常の四角のシミュレーションか酸および重水素化された四角酸共に 600K のら得られた結果であり、b) は重水素化され方が結合長の範囲が広く、結晶全体が高温にた D-SQ のシミュレーション結果から得らなった影響を受けていることが分かる。ただれたものである。

し、水素と重水素との違いは明らかではなく、より詳細な解析が必要である。そこで、本研

また、温度の影響を調べるために図 6 には究では、δH(D)に注目して解析を行った。こ高温時(600K)における酸素間の分布をまとこで、特に重要と考えられる 600K の結果にめた。

ついて議論する。
5

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

これらの物性値の関係性にも注目していく
図 8 には、600K におけるδH(D)の分布をま必要がある。そのため本研究では 2 次元分布とめた。ここで、図 8a は H-SQ を、図 8b はによる解析を実施した。

重水素化された四角酸(D-SQ)のシミュレー

図 9 には、600K における PIMD 計算によっション計算から得られたものである。て得られた二次元分布を示す。このような解析から、水素結合長は D2SQ のほうがより大きa)H-SQく変化する傾向にあることが分かる。また、H/D が中央に分布すると酸素-酸素間距離が短くなることも判明した。また、低温シミュレーションとの比較から、高温(600K)時の D-SQの場合には、この傾向がより強いことを見出した。これらの振る舞いが重水素化による相転移温度の変化に影響を与えていると考えている。
a)

H-SQ

b) D-SQa) D-SQ
図 8 600K のδH(D)ときの 1 次元分布。図中の a)
は H-SQのシミュレーションから得られた結果であり、b) は重水素化された D-SQ のシミュレーシ
ョンから得られたものである。

図 8 から分かるように、600 K は酸素原子間の中央付近に H/D が分布し常誘電体の分布を定性的に再現していると考えられる。また、
D の方が酸素原子間を広く分布しており、質量差により H の方が酸素原子間の中央に分布しやすいことが結果より示唆された。

図 9．600K における H2SQ(左)と D2SQ(右)の水

これまでは 1 次元のみによる解析であった素結合距離(OO)と水素位置(H/D)における二が、現象をより深く理解していくためには、次元分布6

学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

き続き大規模スパコンセンターによる支援

5. 今年度の進捗状況と今後の展望が必要不可欠だと思われる。

今年度の研究進捗については、研究計画に

6. 研究業績一覧（発表予定も含む．投稿中・投稿予定おいて記載した九州大学・情報基盤研究開発センターが所有するスパコン上で動作するは含まない）並列プログラムを整備し、四角酸のシミュレ(1) 学術論文 （査読あり）ーションを実施することに成功した。

2021

四角酸に関する長時間のシミュレーションを多数実施することにより、前節に示した[1] Nuttapon Yodsin, Hiroki Sakagami,

ように、低温時および高温時において、通常

Taro Udagawa, Takayoshi Ishimoto,の四角酸と重水素化された四角酸で異なる

Siriporn Jungsuttiwong, and Masanori

振る舞いをすることを明らかにした。この違Tachikawa,いをより詳細に解析することにより、謎の多い分子結晶系での水素の量子力学的同位体効果現象を明らかにできると考えられる。本"Metal-dopedcarbonnanoconesashighlyefficient

catalystsforhydrogenstorage:

Nuclear quantum effect on hydrogen

研究を実施することによって、そのための手

spillover

法を整備し、研究を遂行するための道筋を構

mechanism",

Molecular

Catalysis, 504, 111486 (2021). DOI:築できたと考えている。
このような研究を実施していくためには、
PIMD のような並列性の高いアルゴリズムをdoi.org/10.1016/j.mcat.2021.111486[2]採用するだけでは十分ではなく、実際のコン

Hiroki

Sakagami,

Masanori

Tachikawa, and Takayoshi Ishimoto,

ピュータ環境で動作することができる並列"Theoretical

性の高いプログラムの開発も同時に行う必studyoftheH/D

isotope effect of CH4/CD4 adsorption

要がある。そのため、研究を遂行することはon a Rh(111) surface using a combined

大変難易度の高いものになる。今回、九州大

plane wave and localized basis sets

学・情報基盤研究開発センターが所有するス

method", RSC Advances, 11, 10253-

パコンを利用する機会が得られ、このような

10257

難易度の高い課題にトライすることが出来(2021).

DOI:

10.1039/D0RA10796D

た。スパコン利用の申請が認められたことは、研究を遂行するために大きな助力となった2020と考えている。
本研究を遂行する上では、温度を変化させ[3] Kazuaki Kuwahata and Masanori

たシミュレーションを行う、多数とプロパテ

Tachikawa,"Nuclear quantum effect onィについて 2 次元分布等を用いて同時に解析

the geometry of NH4+(H2O)", Bull. Chem.

する、といった事柄について良く検討するこSoc. Jap., 93, 1558-1563 (2020).とが重要であった。さらに研究を進めていく

DOI: doi.org/10.1246/bcsj.20200120

ためには、このような多角的な解析が必要不[4] Shumpei Ito, Daisuke Yoshida,

可欠である。今後は、幅広い温度条件や複数

Yukiumi Kita, and Masanori Tachikawa,のプロパティを同時に解析することによっ"First-Principlesて、現象の解明に挑みたいと考えている。引7

QuantumMonte学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

CarloStudies

forDoubleMinimaPrediction

forofInt. J. Quant. Chem., 120, e26275 (5Positronicpages)(2020).

DOI:doi.org/10.1002/qua.26275

Hydrogen Molecular Dianion", J. Chem.
Phys., 153, 224305 (6 pages) (2020).

[9]

DOI: doi.org/10.1063/5.0022673

Yutaro Sugiura, Haruya Suzuki,

TakumaOtomo,

TakaakiMiyazaki,[5] Haruya Suzuki, Takuma Otomo,

Toshiyuki Takayanagi, and Masanori

Ryusei

Tachikawa,

Iida,

YutaroSugiura,"Positron-electron

Toshiyuki Takayanagi, and Masanoricorrelation‐polarization potential

Tachikawa,model"Positron

chloroethenes:
electronbinding

Modeling

inforpositronbindinginpositron-polyatomic molecules", J. Comput.

correlation-polarization

Chem., 41, 1576-1585 (2020). DOI:

potentials

fordoi.org/10.1002/jcc.26200

molecularcalculations", Phys. Rev. A, 102,
DOI:[10] Osamu Kobayashi, Kunihiro Noda,

doi.org/10.1103/PhysRevA.102.052830

Naohiko Ikuma, Dai Shiota, Takayoshi

052830(9pages)(2020).

Ishimoto, and Masanori Tachikawa,[6]

YutaroSugiura,

Toshiyuki"Experimental

Takayanagi,

Masanori

Tachikawa,Analyses of the Oxidation Mechanism"Theoretical calculation of positronof the Poly (arylsilane) Family asannihilationusingthe Side Reaction During the Baking

correlation-Process", J. Phys. Chem. C, 124,

spectrumpositron-electron
polarization

potential",

Int.

J.

andComputational

16149-16158 (2020). (Cover Image)

Quant. Chem., 120, e26376 (5 pages)

DOI:(2020).

doi.org/10.1021/acs.jpcc.0c02416

DOI:

doi.org/10.1002/qua.26376[11] Yusuke Ishida, Haruki Funahashi,[7]Haruki

Tachikawa,

Funahashi,

Masanori

Masanori Tachikawa, and Taro Udagawa,

andUdagawa,"Geometrical H/D Isotope Effect ofTaro"Determining if reaction selectivity

Blue-shifting

Dihydrogen-bonded

can be controlled by the H/D isotope

Clusters", Chem. Lett., 49, 745-748

effect in CH...O interactions", Org.

(2020).

Lett., 22, 9439-9443 (2020). DOI:

doi.org/10.1246/cl.200198

DOI:

doi.org/10.1021/acs.orglett.0c03351
[12]
[8]

Hiroki

Sakagami,

Rathawat

Daengngern,

Osamu

Masanori

Kobayashi, Nawee Kungwan, Chanisorn

Tachikawa, and Takayoshi Ishimoto,

Ngaojampa, and Masanori Tachikawa,"Hydrogen/deuterium adsorption and"Nuclear quantum and H/D isotope

absorption

ineffects on three-centered bonding

palladium using a combined plane

diborane: Path integral molecularwave and localized basis set method",

dynamics simulations", Int. J. Quant.

properties

onand8学際大規模情報基盤共同利用・共同研究拠点 2020 年度共同研究 最終報告書

Chem.,

e26179(2020).

DOI:

doi.org/10.1002/qua.26179[13] Kuniaki Ono, Takayuki Oyamada,
Yukiumi Kita, and Masanori Tachikawa,"Theoretical analysis of the bindingof a positron and pair-annihilationin fluorinated benzene molecules",Eur. Phys. J. D, 74, 1-8 (2020). DOI:
doi.org/10.1140/epjd/e2020-100538-3[14] Yusuke Kanematsu, Hiroyuki SKato, Shinya Yoshimoto, Akira Ueda,
Susumu Yamamoto, Hatsumi Mori, JunYoshinobu,

Iwao

Matsuda,andMasanori Tachikawa, "A computational
examination of the electric-fieldinduced proton transfer along theinterfacehydrogenbond

betweenproton donating and accepting selfassembled monolayers", Chem. Phys.
Lett., 741, 137091 (4 pages) (2020).
DOI:doi.org/10.1016/j.cplett.2020.1
37091[15] Osamu Kobayashi, Tomoki Kato,
TakakoMashiko,

Yohei

Haketa,Maeda,

andMasanori

Hiromitsu

Tachikawa, "Computational simulationof the anion binding association
mechanism contributed to rotation ofpyrrole ring for dipyrrolyldiketone
BF2 complexes", RSC Advances, 10,12013-12024(2020).

DOI:

doi.org/10.1039/C9RA09285D(2) 国際会議プロシーディングス （査読あり）(3) 国際会議発表 （査読なし）(4) 国内会議発表 （査読なし）
9