「先進エネルギー論」/「先進エネルギー特論」
目次
ガイダンス
本講義では、分子シミュレーション、とくに分子動力学法を中心とした理論と 応用について学ぶ。分子シミュレーションがどのように発展してきたかという 歴史的背景をたどりながら、古典力学に基づく運動方程式やポテンシャル関数 の考え方、統計力学との関係など、基礎となる物理化学的概念と計算手法を体 系的に解説する。さらに、材料を主な対象として、実際にPC上で分子動力学シ ミュレーションを実行し、モデル設定から計算実行、結果の出力までの一連の 流れを演習形式で体験する。学生は各自PCを持参し、シミュレーションプログ ラムの基本的な使い方に加えて、Pythonによる科学計算の基礎（数値計算、デー タ処理、可視化など）も身につける。Pythonを用いてシミュレーション結果を 解析し、構造情報、エネルギー、拡散係数や自己相関関数といった動的性質な どを定量的に評価することで、計算データからどのように物理化学的な情報を 引き出し、解釈するかを学ぶことを通じて、分子シミュレーションを材料科学・ 物理化学の研究に活用できる基礎的な素養を養うことを目的とする。
成績評価について
授業への参加度合いとレポートの内容で評価します。
分子動力学計算の理論と実行
分子動力学計算の理論
分子動力学計算の基礎理論と歴史、LAMMPSコードについてを説明します。
自前のPCでLAMMPSによるMD計算を実行し、ログの確認とMD計算の結果をVMDで 可視化できるようになることを目標とします。
LAMMPS
Windows
LAMMPSはUSのSandia National Laboratoryで無料で配布されている汎用分子 動力学計算プログラムです。LAMMPSの実行ファイルを入手し、コマンドラ インから計算方法を記述した入力ファイルを指定することで、MD計算を 行うことができます。MD計算のための環境を整えるのがやや面倒なので、簡 単に実行できるよう実行環境を大窪が作成しています。以下の手順に従って 計算環境を作ります。
- LAMMPSのWindows用バイナリファイル(lammps_win64.zip)をダウンロード して展開する。
展開した
lammpsフォルダをC:\lammpsに配置する。
Windowsのスタートメニューを開いて「command」を検索しコマンドプロン プトを実行します。
以後全てコマンドプロンプト上でコマンドラインにて作業を行います。最低 限、次のコマンドを覚えておけば大丈夫です。(矢印の↑でコマンドの履歴 を呼び出せます)
cd folder_name folder_nameに移動 cd .. 1つ上のフォルダに移動 dir folderに存在するファイル一覧を表示 (macOSではdirの代わりにls) Ctrl-c Ctrl-cで実行中のプログラムを強制終了
次のコマンドで C:\lammps に移動し、サンプル計算を実行します。
cd /d C:\lammps lmp -in sample.in
次のように計算結果が表示され、 ERROR が出ずに計算が終了し、
作業フォルダに log.lammps が生成されることを確認してください。
細かい数字は上の出力結果と異なっていても問題ないです。
macOS
- lammps_osx_arm64.zip (Apple Silicon: arm64)
- lammps_osx_x86.zip (Intel: x86_64)
- ZIPを展開した後、ターミナルで展開先のフォルダを開きます。
次のコマンドで
lammpsフォルダに移動し、サンプル計算を実行します。 実行ファイルのパスを示す./を忘れないようにしてください。 Windowsの場合と同様に、ERRORが出ずに計算が終了し、log.lammpsが生成されることを確認してください。
cd lammps ./lmp -in sample.in
Linux
compile済の実行ファイルが配布されています。 https://docs.lammps.org/Install_linux.html
VMD
VMDはLAMMPSから得られた計算結果を可視化するプログラムです。公式サイトでユーザー登録してから自分のOSに合ったファイルを入手してください。
- Windows: 展開した
VMDフォルダ内のvmd.exeを実行します。 日本語を含まないパスに配置してください。 - macOS: 展開したVMDアプリ（=.app=）を「アプリケーション」フォルダなどに 移し、起動します。Apple Silicon用とIntel用を選び分けてください。
- Windows用: vmd_win32.zip
- macOS用 (Apple Silicon: arm64): vmd_osx_arm64.tgz
- macOS用 (Intel: x86_64): vmd_osx_x86.zip
Pythonによるデータ解析
Pythonを用いて、MDから出力される統計データや座標データを解析できるようになることを目標とします。
Python実習
Google Colabを利用してPythonの基礎とデータ解析を説明します。各自、大学から配布されているGoogleアカウントで、ログインしてGoogle Colabを利用してください。
Web版の講義資料も参照してください。
液体ArのMD
論文 Correlations in the Motion of Atoms in Liquid Argon.pdf
- Author: A. Rahman (Argoone National Lab.)
- Physical Review A, 136, 405-411 (1964)
- Abstruct: A system of 864 particles interacting with a I.ennard-Jones potential and obeying classical equations of motion has been studied on a digital computer (CDC 3600) to simulate molecular dynamics in liquid argon at 94.4 K and a density of 1.374 g/cm3. The pair-correlation function and the constant of self-di6usion are found to agree well with experiment; the latter is 15%lower than the experimental value. The spectrum of the velocity autocorrelation function shows a broad maximum in the frequency region of the Van Hove function G,(r,t) attains a maximum departure from a Gaussian at about s=3 X10 "sec and becomes a Gaussian again at about 10 "sec. The Van Hove function G(r, t) has been compared with the convolution approximation of Vineyard, showing that this approximation gives a too rapid decay of Gd, (r,t) with time. A delayed-convolution approximation has been suggested which gives a better 6t with Gz(r, I); this delayed convolution makes Gz(r, t} decay as t4 at short times and as I at long times.
液体Arの分子動力学計算
インプットファイル
- 密度は1.374 g/cm3(論文と同じ)とします。入力ファイル: Ar.zip
- インプットファイル理解します。インプットファイルの解説: Ar_in_data.pdf
- lammpsを実行します。-inでインプットファイル、-logで出力の保存ファイルを指定します。
- 標準出力にエネルギーや温度等の系の情報が出力されます(Ar.logに同じ内容が保存されている)。Ar.inを編集して温度や計算時間等を変えてMD計算を実行します。
- 動径分布関数と平均自乗変位を計算する入力ファイル: Ar_nvt.zip
アルカリハロゲン溶融塩のMD
講義資料
参考論文
Interatomic Distances in Crystals of the Alkali Halides.pdf
- J. Chem. Phys., 1, 643-646 (1933)
The constants of the empirical repulsion potential assumed by Born and Mayer between two ions, which may be termed the ionic "radii," are calculated for the alkali and halide ions from the lattice constants of the alkali halides. By using 9 such radii the 20 lattice constants can be recalculated nearly within the probable experimental error of their determination. The good agreement is regarded as support for the assumptions involved in the calculations.
Ionic sizes and born repulsive parameters in the NaCl-type alkali halides—I.pdf
- Author: F. G. Fumi and M. P. Tosi
- J. Phye. Chem. Solids, 25, 31-43 (1964)
- Abstruct: The Born model of ionic solids is shown to permit the complete determination from solid-state data of crystal radii for the alkali and halogen ions in the individual alkali halides with the NaCl structure, which are directly comparable with the experimental crystal radii deducible from the X-ray maps for the electron distribution in the crystal.
Ionic sizes and born repulsive parameters in the NaCl-type alkali halides—II.pdf
- Author: F. G. Fumi and M. P. Tosi
- J. Phys. Chem. Solids, 25, 45-62 (1964)
- Abstruct: procedure described in I for the determination of the crystal radii of the ions in the individual NaCl-type alkali halides from solid-state data by the Born model is applied, adopting a generalized Huggins-Mayer form for the Born repulsive energy (i.e. allowing the hardness para- meter to vary from salt to salt).
Calculation of the melting point of NaCl by molecular simulation.pdf
- Author: J Anwar and D. Frenkel
- J. Chem. Phys., 118, 728-735 (2003)
- Abstruct: We report a numerical calculation of the melting point of NaCl. The solid–liquid transition was located by determining the point where the chemical potentials of the solid and liquid phases intersect.
イオン液体のMD
講義資料
イオン液体のMD計算の論文(カチオン)とアニオンを参考にポテンシャル関数等 を解説した後、1-ethyl-3-methylimidazolium (c2c1im)と trifluoromethylsufate (ntf2)からなるイオン液体のMD計算を行います。構造 やダイナミクスを評価するとともに粘性係数を計算します。
実習: イオン液体のMD
分子系のインプットファイルは、各分子の結合距離や角度、二面角等を指定す る必要があります。分子系のインプットファイルを手作業で作成することは、 あまりにつらいので配布したスクリプトで作成します。c2c1imとntf2が27個ず つからなるセルのインプットファイルを以下のコマンドで生成します。
make_il.py --c2c1im 27 --ntf2 27 -o emim-ntf2
このコマンドで以下の2つのファイルが生成されます。
- emim-ntf2.data
- emim-ntf2.in
emim-ntf2.dataは、原子が重ならないように大きめの単位胞を初期構造として います。短いステップで平衡構造を得るために2,000 KのNVT計算でアニーリン グし、deformを使って1.5 g/cm3まで体積を減少させならがら300Kまで温度 を下げます。その後、1 atmで300 KのNPT計算を行って平衡構造を得ます。