分子動力学計算の実践講習 −溶融塩とイオン液体−

目次


日時: 2019年11月19日(火) 9:00〜18:00 【事前講習会11月18日(月)13:00〜18:00】
場所: 千葉大学西千葉キャンパス (東京駅から最寄りのJR西千葉駅まで1時間弱)
会場: 千葉大学ベンチャービジネスラボラトリ2階会議室 (JR西千葉駅から徒歩5分)
担当: 千葉大学大学院工学研究院 大窪 貴洋

申し込み方法や参加費の支払い方法などは、(社)電気化化学会溶融塩委員会 溶融塩化学講習会のページを御覧ください。

注意事項

11/19(火)の講習会は、持ち込みPCの LinuxまたMacOSX上でLAMMPSを実行します。可視化はVMDとVESTAで行います。 可視化のためのVMDやVESTA、解析のためのPythonを本講習会の前にインストー ルしておいてください。一部、解析のためにPythonを使用します。Linux環境 のセットアップやLAMMPSのインストール等は取り扱いませんので、自力での環 境構築が難しい場合は前日11/18(月)の事前講習会にPCを持参して参加してく ださい。事前講習会ではMacOSXまたはWindows10上のWSLでLinux環境を構築し てLAMMPSが実行できるようセットアップします。またUNIXのコマンドライン操 作に不安のある人も事前講習会に参加してコマンドライン操作でのLAMMPSの実 行に慣れておいてください。

はじめに

材料の分子構造や物性を予測するため、計算機シミュレーションがあたりまえ のように使われてきています。特に溶融塩やイオン液体のように高温下や不活 性雰囲気での取扱いを必要とし、実験することが容易でない材料の特性や分子 構造を調べるために計算機シミュレーションは有用です。溶融塩分野での計算 機シミュレーションは、原子・電子構造解析から反応炉の設計まで幅広く用い られていますが、本公演では溶融塩の物性予測や分子構造の解析に有用な分子 動力学計算を取り扱います。本講習会は、分子動力学計算の初学者を対象とし、 参加者が次の項目を達成できるよう講習します。

  • 分子動力学計算の基本原理を理解する。
  • 論文や書籍の分子動力学計算を自分でトレースできるようになる。
  • 組成や温度をパラメータとしたシミュレーションを行い溶融塩の物性を予測する。

現在、分子動力学計算を行うソフトウェアは無料のものから有料のものまで多 数存在します。シミュレーションの実行環境もデスクトップPCから大型の並列 計算機からまで多岐にわたります。本講習会の達成目標を考慮すると、試行錯 誤を何度も繰り返してシミュレーションを行い、理解を深めてもらった方が望 ましいので、自前のPCを使って講習することとします。また分子シミュレーショ ンコードとして、学術研究や企業の研究開発でもよく用いられているLAMMPSを 用います。LAMMPSは、デスクトップPCから大型の並列計算機まで環境を選ばず 実行できますので、自前のPCで慣れてきたら大型の計算機上へ実行環境を移行 することもスムーズかと思います。またGUIを搭載した以下の分子動力学計算 用プリ・ポスト処理ソフトウェアは分子動力学計算としてLAMMPSを利用して います。

このようなプリ・ポスト処理ソフトウェアを利用して計算を行う場合でも LAMMPSをよく理解してシミュレーションを行うことは、動作原理を理解する上 でも有用と思います。またプリポスト用ソフトウェアで対応していないLAMMPS の機能や、最新のLAMMPSバージョンを利用したい場合は、LAMMPSを直接実行す る必要がありますので、LAMMPSの使い方を覚えておいて損はありません。

本講習会では、LAMMPSのインプット・アウトプットファイルの形式を理解し、自前のPCで結晶構造や分子構造、ガラス組成からMD計算をスタートする方法を説明します。また大型の並列計算機が一般にLinux等のUNIX環境であることから、Linux環境でのシミュレーションを講習します。講習会では、LinuxやLAMMPSのセットアップは行いませんので、自力での構築が難しい場合は、事前講習会に参加してください。

講習会スケジュール

進行具合を確認しながら進めますので、タイムテーブルはあくまで予定です。

9:00 9:05 講習会はじまりの挨拶
9:00 10:00 分子動力学計算、LAMMPSコードについて (講習)
10:00 11:00 LAMMPSの実行、インプット・アウトプットについて (講習・実習)
11:00 12:00 液体Arの分子動力学計算 (講習・実習)
12:00 13:00 昼休み
13:00 14:30 溶融塩のシミュレーション (講習・実習)
14:30 16:00 イオン液体のシミュレーション (講習・実習)
16:00 17:30 インプットファイルの作成、シミュレーションの実行 (自由課題)

分子動力学計算の基本

分子動力学計算の基礎的事項の講義

  • 分子動力学計算の歴史
  • 速度Verlet法
  • 周期境界条件
  • ポテンシャル関数
  • アンサンブル
  • 自己相関関数

LAMMPSコードについて

  • コードの特徴
  • ソースコードの構成

液体Arの分子動力学計算

液体Arの分子動力学計算を行った論文を取り上げ内容を解説した後、論文と同 様のMD計算を行い、液体Ar構造の可視化、密度評価、動径分布関数による構造 評価、ダイナミクスとして平均自乗変位から自己拡散係数を求めます。

初期構造の作成

Arの密度からセルを作る。ポテンシャル関数の変数も定義する。

Liquid Ar (rho=1.374g/cm3), Phys. Rev. 136, A405 (1964)

1 atoms

1 atom types

0 3.6412664191749493 xlo xhi
0 3.6412664191749493 ylo yhi
0 3.6412664191749493 zlo zhi

Masses

1 39.948


Pair Coeffs # lj/cut

1 0.23846475655872523 3.4

Atoms

1   1     0.0 0.0 0.0

平衡計算と動径分布関数と平均自乗変位の計算

  • 10x10x10のセル(1000個のAr)としNVT(94.4K)とする。
  • 3 psの計算で平衡状態を得る。
  • 12 Åまでの動径分布関数を計算する。
  • 100 fs毎に平均自乗変位を計算する。
units          real
pair_style     lj/cut 7.65
read_data      Ar.data
replicate      10 10 10
thermo         10
thermo_style   custom step temp etotal press vol density
velocity       all create 94.4 777
fix            1 all nvt temp 94.4 94.4 10
dump           1 all custom 1 Ar_nvt.lammpstrj id element x y z
dump_modify    1 element Ar sort id
timestep       10
run            300

comm_modify    cutoff 14
compute        a all rdf 240 * * cutoff 12
fix            2 all ave/time 1 300 300 c_a[*] file Ar_nvt.gr mode vector
compute        b all msd
fix            3 all ave/time 1 1 10 c_b[*] file Ar_nvt.msd title1 "#DT=10fs"
run            300
  • 平衡状態の確認
lmp_logplot.py Ar_nvt.log -k TotEng Temp -w 30 60 -o log.png

log.png

図1: エネルギーと温度の変化

snapshot.png

図2: 初期構造と6.0 ps後のスナップショット

構造とダイナミクス評価

  • 動径分布関数と積算配位数
lmp_grplot.py Ar_nvt.gr -o gr.png

gr.png

図3: Arの動径分布関数(94.4 K)

  • 平均自乗変位と拡散係数
lmp_msd.py -w 1 3 Ar_nvt.msd -o msd.png

msd.png

図4: Arの平均自乗変位(94.4 K)

練習

  • 1atmで50K, 100K, 150K, 200K, 250KのNPT計算を行う。
  • 各温度での平衡になるステップ数を確認したら、平衡後のMD計算に基づいて密度および動径分布関数と拡散係数を評価する。

固体NaClと溶融NaClの分子動力学計算

MD計算による溶融NaClの熱伝導率を計算した論文を取り上げポテンシャル関数等を解説した後、NaClのMD計算を行います。構造やダイナミクスを評価するとともに熱伝導率を計算します。

初期構造の作成

NaClの結晶構造データ(Crystallographic Information File; CIF)から初期 構造を作成します。

COMMENT

8 atoms
2 atom types

0.0  5.620  xlo xhi
0.0  5.620  ylo yhi
0.0  5.620  zlo zhi

PairIJ Coeffs

1 1 6.081 0.31700 2.340 24.18024 -11.5120
1 2 4.864 0.31700 2.755 161.2016 -200.02100
2 2 3.649 0.31700 3.170 1669.588 -3352.8700

Masses

1  22.989768
2  35.4527

Atoms

1   1   +1.0    0.000000    0.000000    0.000000  # Na
2   1   +1.0    0.000000    2.810000    2.810000  # Na
3   1   +1.0    2.810000    0.000000    2.810000  # Na
4   1   +1.0    2.810000    2.810000    0.000000  # Na
5   2   -1.0    2.810000    2.810000    2.810000  # Cl
6   2   -1.0    2.810000    0.000000    0.000000  # Cl
7   2   -1.0    0.000000    2.810000    0.000000  # Cl
8   2   -1.0    0.000000    0.000000    2.810000  # Cl

初期構造から次のステップで計算を行って、300Kと1200KのMD計算を行います。

variable       RUN equal  1000
variable       DT equal 10.0
variable       cutoff equal 10            
variable       nbin equal ${cutoff}/0.01
variable       title format DT "#DT=%.8gfs;"
variable       elems string "#ELEMENTS=Na-Cl;"

units          real
atom_style     charge
pair_style     born/coul/long ${cutoff}
read_data      NaCl.data 
replicate      4 4 4
group          Na type 1
group          Cl type 2
kspace_style   ewald 1.0e-4
thermo	       100
thermo_style  custom step temp etotal press vol density
timestep       ${DT}

# 300 K
variable        T equal 300
velocity        all create ${T} 777
fix	        1 all npt temp ${T} ${T} 10 iso 1 1 100
dump            1 all custom 10 NaCl.lammpstrj id type element x y z
dump_modify     1 element Na Cl sort id
run             ${RUN}
compute         1 Na msd 
compute         2 Cl msd 
compute         3 all rdf ${nbin} 1 1 1 2 2 2
fix             2 all ave/time 1 1 5 c_1[*] file Na_${T}K.msd title1 ${title}
fix             3 all ave/time 1 1 5 c_2[*] file Cl_${T}K.msd title1 ${title}
fix             4 all ave/time 1 ${RUN} ${RUN} c_3[*] file NaCl_${T}K.gr mode vector title1 ${elems}
run             ${RUN}
uncompute       1
uncompute       2
uncompute       3
unfix           2
unfix           3
unfix           4

# 2000K->1200K
variable        T equal 1200
fix	        1 all npt temp 2000 2000 10 iso 1 1 100
run             ${RUN}
fix	        1 all npt temp 2000 ${T} 10 iso 1 1 100
run             ${RUN}
# 1200K
fix	        1 all npt temp ${T} ${T} 10 iso 1 1 100
run             ${RUN}
compute         1 Na msd 
compute         2 Cl msd 
compute         3 all rdf ${nbin} 1 1 1 2 2 2
fix             2 all ave/time 1 1 5 c_1[*] file Na_${T}K.msd title1 ${title}
fix             3 all ave/time 1 1 5 c_2[*] file Cl_${T}K.msd title1 ${title}
fix             4 all ave/time 1 ${RUN} ${RUN} c_3[*] file NaCl_${T}K.gr mode vector title1 ${elems}
run             ${RUN}
write_data      NaCl_1200K.data pair ij

MD計算でのエネルギーと温度、密度変化を確認します。

lmp_logplot.py NaCl.log -k TotEng Temp Density

log.png

図5: エネルギー、温度、密度変化

snapshot.png

図6: 初期構造から温度を変えて計算して得られた構造のスナップショット

構造とダイナミクス評価

動径分布関数をプロットします。

lmp_grplot.py  NaCl_300K.gr -o NaCl_300K.gr
lmp_grplot.py  NaCl_1200K.gr -o NaCl_1200K.gr

gr.png

図7: NaClの動径分布関数, 300K(上)と1200K(下)

平均自乗変位をプロットして、拡散係数と温度因子を計算します。

lmp_msd.py Na_1200K.msd Cl_1200K.msd -w 5 10
lmp_msd.py Na_300K.msd Cl_300K.msd -w 5 10 --order 0

msd.png

図8: NaClの平均自乗変位, 300K(上)と1200K(下)の

熱伝導率の計算

保存した1,200Kでのデータを読んでフラックスの計算します。

variable        T equal 1200       # temperature
variable        DT equal 2         # timestep
variable        V  equal vol       # volume
variable        RUN equal 10000    # run period
variable        OUTFLUX string NaCl_${T}K.flux    

units           real
pair_style      born/coul/long 10
kspace_style    ewald 1.0e-4
atom_style      charge
read_data       NaCl_${T}K.data

variable        DT_ format DT "DT=%.8gfs"
variable        T_ format T   "T=%.8gK"
variable        V_ format V   "V=%.8gAng3"
variable        title string  "#${DT_};${T_};${V_};[Kappa_calc.]"

timestep        ${DT}
velocity        all scale ${T}
fix             1 all nve
# compute energy flux
compute         myKE all ke/atom
compute         myPE all pe/atom
compute         myStress all stress/atom NULL virial
compute         flux all heat/flux myKE myPE myStress
variable        Jx equal c_flux[1]/$V
variable        Jy equal c_flux[2]/$V
variable        Jz equal c_flux[3]/$V
# save flux data every step
fix             2 all print 1 "$(step) ${Jx} ${Jy} ${Jz}" &
                title ${title} screen no file ${OUTFLUX}
thermo	        100
thermo_style    custom step temp etotal press vol density
run             ${RUN}

0.3 psの時間windowで相関関数を計算して積分し熱伝導率を求めます。

lmp_corr.py NaCl_1200K.flux  -s 1 -w 0.3

corr.png

図9: 相関関数と相関関数の積分から求めた熱伝導率

イオン液体の分子動力学計算(予定)

イオン液体のMD計算の論文(カチオン)アニオンを参考にポテンシャル関数等を解説した後、1-alkyl-3-methylimidazoliumとtrifluoromethylsufateからなるイオン液体のMD計算を行い、構造やダイナミクスを評価するとともに粘性係数を計算します。

初期構造の作成

平衡計算

構造とダイナミクス評価

粘性係数の計算

自由課題の分子動力学計算

参加者の計算したい材料(ガラス、結晶、イオン液体等)のインプットファイルを作成し、計算行います。時間の許す限り担当者がアドバイスしますので何でもお尋ねください。

事前講習会

MacOSXまたはMay 2019以降のwindows10にupdateしたPCを持参してください。windowsメニューの設定からバージョン情報を確認してバージョンが1903以上だったらOKです。windows10のupdateは時間がかかるので、必ず事前講習会に参加するまでにupdateしておいてください。

windws_version.png

LAMMPSをセットアップして、自前のPCで実行する環境を構築します。またプリ /ポスト処理でPythonによるデータ処理を行うため、Pythonのセットアップも 行います。エディタは個人の好きなものを選択して良いですが、UNIXの改行コー ドでUTF-8でエンコードしたファイルを編集できるものを使ってください。特 に何もないようならEmacsをおすすめします(EmacsのLammpsモードを使えば、 インプットファイルが編集しやすいため)。

スケジュール

13:00 14:00 Linux, LAMMPS, Pythonのセットアップ
14:00 15:00 UNIXコマンド、簡単なシェルスクリプトの練習
15:00 16:00 エディタ(Emacs)の使い方
16:00 17:00 LAMMPSの実行

WSLのインストール

  • Windowsのスタートメニューからコントロールパネルを検索して実行します。

    control_panel.png

  • プログラムと機能からWindowsの機能の有効化または無効化をクリックします。

    windows_add.png

  • Windows Subsystem for Linuxにチェックを入れてOKを押します。これで Windows10上でLinuxを動作させることができるようになります。

    windows_func.png

Debianのインストール

  • windows storeをスタートメニューから選びます。

    windows_store.png

  • Debianを検索してインストールします。これでLinuxのディストリビューションの1つであるdebianがインストールされます。

    debian.png

VMDのインストール

計算結果を可視化するためにVMDをインストールします。以下のURLに接続し てWindows用のファイルをダウンロードします。最初にダウンロードする時に、 Usernameの登録が求められるので指示に従います。 https://www.ks.uiuc.edu/Research/vmd/

VESTAのインストール

結晶構造の共通フォーマットであるCIFファイルを可視化したり対称性を操作するためにVESTAをインストールします。 https://jp-minerals.org/vesta/jp/

Pythonのインストール

データ処理用にPythonをセットアップします。インストール用シェルスクリプトをダウンロードして実行します。

curl -O https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh

インストール先を聞いてきますが、デフォルト($HOME/miniconda3)で問題ありません。 Python本体をインストールしたら必要なモジュールをインストールします。

conda install numpy matplotlib

その他

Lammpsのコンパル

標準でインストールされていない機能を使いたい、Lammpsを改造して使いたいというときは、Lammpsのソースコードを取得して自前でコンパルする必要があります。 Lammpsのコンパルとリンクに必要なコンパイラ(g++)と高速フーリエ変換ライブラリ(fftw3)をaptでインストールしてgitでLammpsのソースコードを取得します。intel製のコンパイラやMKLライブラリを購入している人は、そちらでコンパイル・リンクした方が早いです。

PCをまたいだ複数のノードで計算した場合は、mpiを導入する必要がありますが、こちらの構築は割愛します。研究室で構築した並列計算機のメモが少しは役に立つかもしれません。

apt install g++ fftw3-dev git

gitを使ってlammpsのソースコードを取得します。

git clone -b stable https://github.com/lammps/lammps.git lammps

ソースコードを取得したら、lammps/srcに移動します。

cd lammps/src

必要な拡張機能のソースコードを追加します。

make yes-RIGID
make yes-USER-REAXC
make yes-MISC
make yes-KSPACE
make yes-MOLECULE

MAKEファイル以下にあるコンパイル用のMakefileを編集します。ewald法等の長距離クーロン相互作用を計算するために、フーリエ変換ライブラリ(fftw3)を利用するので、ライブラリの場所をMakefile中で指定します。複数コアを使った計算を行うためにopenmp付きでコンパイルします。openmpiを使わない場合は、dummy用のライブラリを作っておく必要がありまうす。STUBSに移動してライブラリを作成しておきます。

cd STUBS
make CC=g++-mp-8
cd ..
make mac_serial -j8

コンパイルします。Makefile.hogeならばlmp_hogeの実行ファイルが生成されます。-jでコンパイルの並列処理数を指定できます。

本講習会について

本講習会は電気化学会溶融塩委員会の主催で開催されます。 http://msc.electrochem.jp/index.html

著者: 大窪 貴洋

Created: 2019-07-15 月 00:13

Validate