どぼるざーくとゆく、高分子シミュレーション

lammpsやOCTAなどの高分子シミュレーターについて四苦八苦しながら学んでいくブログです。Twitterアカウント@lammps_octa

【お遊び】水分子と愉快なモデル達

f:id:toal_peg_simu:20200121215809p:plain

MD戦隊、水レンジャー!!!

はじめに

みなさんこんにちは!

前回の記事が思ったよりも長くなってしまったことと、内容が堅苦し過ぎたので今回の記事はネタに走りたいと思います(笑)

テーマは水分子です!!私たちの周りにある『水』ですが、シミュレーション界隈ではどのような扱いを受けているのでしょうか??早速見てみましょう!

 

 

 

いでよ!水分子のモデル

シミュレーション界隈に言ってみた

今やっている高分子の研究は溶媒中での挙動が重要になってくるので、シミュレーション界隈の人に「水分子を扱ってるんですよ~」と話すんですが、ほとんどの人が

「水分子は鬼門、触れてはいけない」

「学部生で水分子を扱うのは大変では?」

「全然安定しないよね~」

という辛辣なコメントをもらいます(苦笑)

実は水分子は構造がシンプルですが水素と酸素の間で量子効果が働くため、ポテンシャルを当てはめるだけでは水分子を表現することはできません(=゚ω゚)ノ

そこで水分子をMDシミュレーションで用いる場合には、用いる系に合うようなモデルを選択して使うことになります。本当にたくさんの種類があるので、今回は4つ紹介します!

ちなみに参考にした論文は『Flexible simple point-charge water model with improved liquid-state properties』というものです。

(https://aip.scitation.org/doi/full/10.1063/1.2136877?ver=pdfcov)

 

SPC

f:id:toal_peg_simu:20200121221634p:plain

このモデルは各粒子に電荷があり結合長・結合角が変化しないモデルです。一番簡単なモデルではありますが、定常状態(300K,1atm)では水の物性を比較的高精度に再現することが出来ます!

一般的に水分子は生体たんぱく質などのMD計算で必要になるのですが、相互作用の要素が1つ増えるごとに解析時間が最大50%も増加することが知られています。そのため出来るだけ少ないポテンシャルで水分子を表現できた方がいいんですよね。

 

TIP3P

f:id:toal_peg_simu:20200121221909p:plain

こいつもSPCモデルと同様に各粒子に電荷があり結合長・結合角が変化しないモデルです。SPCモデルとは電荷と結合長のパラメータが微妙に異なります。

ボンドの結合長 SPC:1.0Å TIP3P:0.9527Å

うん、、、、、めっちゃ少しやん!!もはや10のマイナス12乗レベルでしか誤差ないやん( ̄д ̄) でもこの小さな差で物性値に大きな差が出るんですよね。例えば拡散係数はSPCよりもTIP3Pの方が1.3倍も大きくなってしまいます。

ここで察して頂けたと思いますが水分子めっちゃセンシティブ。恋する乙女並みにパラメータに敏感なのです♡つまり水分子を実際の分子に近づけるには、かなり桁数を増やして計算をする必要がありそうですね(白目)

ちなみにこのモデルは生体分子のシミュレーションによく使われているらしいです。その理由としてはCHARMMと呼ばれる分子動力学シミュレーターの力場がTIP3Pをもとに作成されたからだそうです。解説記事があるので、気になる人はそちらをチェック↓

masa-cbl.hatenadiary.jp

SPC/Fd

f:id:toal_peg_simu:20200121221919p:plain

このモデルはSPCモデルを参考に作られたモデルですが、結合長・結合角が変化するという点が今までのモデルと異なります!実際の水分子内の結合も変化しているので、より現実に即したモデルと言えるのでしょう。

ただ現実と同じように結合長を可変にしたからといって物性値が水に近づくわけではないんですね。。。なんなら結合長を変えないモデルの方が実験に近い物性値もあるので計算時間を無駄に増やしただけなのではと思ってしまいます(シミュレーション難しい)

SPC/Fw

f:id:toal_peg_simu:20200121221944p:plain

最後に紹介するのはこちら!SPC/Fwモデルです。

え、SPC/Fdモデルと変わらなくない?と思った人もいるかもしれませんが、少しは変

えてるんです!!!

結合長(Å):1.0→1.012

結合角(θ):109.5→113.24

 いや、ほぼ誤差じゃん??と最初見た時には私も思いました(笑)

ただこの少しのパラメータの変更だけで、拡散係数・水素結合継続時間・粘性などの物性値が現実のモデルの物性値に近づいたそうです。特に拡散係数や粘性は生体分子のシミュレーションに大きな影響を与えるため、このモデルを使った方が生体で起きる反応をより正確に記述できるのではないでしょうか。

ちなみにこのモデルは私がシミュレーションしたことがあるので動画を挙げておきますね!


【お遊び】水分子と愉快なモデル達~SPC/Fwモデル~

さいごに

いかがでしたか?

このように身近な『水』をシミュレーションに導入するといっても、モデルの選択や計算負荷の増大、パラメータの敏感さなど考慮する点が多く存在するということが分かってもらえたでしょうか?

ちなみに参考にした論文を読むことが出来れば、lammpsでも簡単に再現することができるので興味を持った人は是非やってみてください!

【lammps】pythonで高分子を書いてみよう!!

はじめに

皆さん、こんにちは!

明日試験なんですけど、現実逃避したいのでpythonを使った分子の初期構造の書き方を解説していきたいと思います!

lammpsを使っている人は既に他の言語をやっている人も多いと思いますが、python使えるようになると出力データの解析も簡単に行えるようになるので、ぜひチャレンジしてみましょう!

 

pythonとは?

pythonとは機械学習などが初心者でも簡単に出来るプログラミング言語の一種です。pythonJapanのHPでは以下のように解説されています。

とてもクリーンで読みやすい文法

強力な内省(イントロスペクション)機能

直感的なオブジェクト指向  ・・・・・・・

(https://www.python.jp/pages/about.html)

 私がlammpsで自動化を行う時にはPerlpythonを利用しているのですが、その理由としては、「直感的に書ける」「多くの人が使ってるので、インターネットや友人に疑問点を聞きやすい」「ライブラリーが多くファイル操作や数値解析がしやすい」などがあります。

今回は初期構造を作成するプログラムについて解説しますが、今後の記事で出力データの解析やファイル操作についても解説していきたいと思います!

(※pythonを使う環境は特に指定しませんが、『math』『random』パッケージが最低でも使えるようにしておいてください。)

 

 

pythonのコードを解説

今回のプログラムで得られる初期構造

今回のプログラムでは以下のような初期構造を作成することができます。

・重合度nの粗視化高分子鎖がセルに1本ある

 

ここで言う粗視化高分子鎖とは、側鎖が無く粒子の種類が1種の高分子のことです。重合度nは各自で書き換えることができるので、色々試してみてください!

コード全体

#使用するパッケージ
import math
import random


#各種変数設定-----------------------------------------------
#原子量初期設定
mass_H=1.00
mass_C=12.0
mass_O=16.0

#粒子の重さ
mass_g=mass_C*2+mass_H*4+mass_O

#重合度からセルサイズを決定--------------------------------
print('重合度をどうしますか?(nは4以上)')
#高分子鎖の重合度
n=int(input())
c_leng=(n-1)*3.30
cut_leng=12

#セルサイズ
print('この場合の高分子の最大長さは'+str(c_leng)+'です。セルの一辺の長さはどうしますか?')
print('※今回のカットオフの長さは'+str(cut_leng)+'です。')
cellsize=int(input())
x=cellsize
y=cellsize
z=cellsize


#各種設定-------------------------------------------------------
#高分子の色んな数
c_bond_n=n-1
c_angle_n=n-2
c_dihed_n=n-3

atoms=n
atom_types=1
bonds=c_bond_n
bond_types=1
angles=c_angle_n
angle_types=1
dihedrals=c_dihed_n
dihedral_types=1



print('#各種結合の設定')
print(str(atoms)+' atoms')
print(str(atom_types)+' atom types')
print(str(bonds)+' bonds')
print(str(bond_types)+' bond types')
print(str(angles)+' angles')
print(str(angle_types)+' angle types')
print(str(dihedrals)+' dihedrals')
print(str(dihedral_types)+' dihedral types\n')



#初期セルサイズ-------------------------------------------------------
print('0 '+str(x)+' xlo xhi')
print('0 '+str(y)+' ylo yhi')
print('0 '+str(z)+' zlo zhi')
print(' \n')


#質量設定-------------------------------------------------------
print('Masses\n')
#高分子鎖の原子設定
print('1 '+str(mass_g))

print(' \n')


#Atoms-------------------------------------------------------
print('Atoms\n')

#変数
#Noの為の変数
no=1
#炭素鎖の奥行方向の座標
c_x0=x/2
#炭素鎖の軸方向の座標
c_y0=y/10
#炭素鎖の高さ方向の座標
c_z0=z*3/4
#粒子のの配置間隔
#wantleng=
#step=n*3.3/wantleng
step=1.0


#高分子鎖の原子位置
for i in range(n):
    #左の炭素
    print(str(no)+' 1 1'+' '+str(c_x0)+' '+str(c_y0)+' '+str(c_z0)+' 0 0 0')
    no+=1
    c_y0+=step

print(' \n')

#Bonds-------------------------------------------------------
print('Bonds\n')

#Bondの変数
b_no=0
#高分子鎖のBondをつなぐための変数
b_n1=1
b_n2=2


#繰り返し部分
for i in range(c_bond_n):
    b_no+=1
    print(str(b_no)+' 1 '+str(b_n1)+' '+str(b_n2))
    b_n1+=1
    b_n2+=1
print(' \n')


#Angles----------------------------------------------------------------------------------
print('Angles\n')

#Bondの変数
a_no=0
#高分子鎖のangleをつなぐための変数
a_n1=1
a_n2=2
a_n3=3

#繰り返し部分
for i in range(c_angle_n):
    a_no+=1
    print(str(a_no)+' 1 '+str(a_n1)+' '+str(a_n2)+' '+str(a_n3))
    a_n1+=1
    a_n2+=1
    a_n3+=1

print(' \n')


#Dihedrals------------------------------------------------------------------------------------
print('Dihedrals\n')

#dihedralの変数
d_no=0

#高分子鎖のdihedralをつなぐための変数
d_n1=1
d_n2=2
d_n3=3
d_n4=4

#繰り返し部分
for i in range(c_dihed_n):
    d_no+=1
    print(str(d_no)+' 1 '+str(d_n1)+' '+str(d_n2)+' '+str(d_n3)+' '+str(d_n4))
    d_n1+=1
    d_n2+=1
    d_n3+=1
    d_n4+=1

 コードの解説

それではこのコードの重要ポイントを解説したいと思います。

#使用するパッケージ
import math
import random

 ここの部分ではプログラムの中で使うパッケージを『import ほげほげ』で指定しています。もし数値計算などを行う場合にはpandasなどを入れる必要があります。

#各種変数設定-----------------------------------------------
#原子量初期設定
mass_H=1.00
mass_C=12.0
mass_O=16.0

#粒子の重さ
mass_g=mass_C*2+mass_H*4+mass_O

 ここでは粗視化粒子の重さを規定しています。

『#原子量初期設定』の項目では粗視化粒子内に存在する原子の原子量を指定しています。今回の粗視化高分子はポリエチレングリコールを想定しているので、水素・炭素・酸素の変数を作っています。(※もし他の原子が入る粗視化粒子をする場合は自分で書き加えましょう)

『#粒子の重さ』の項目では今回のモデルが『-(CH2CH2O)n-』を1つの粗視化粒子を想定しているため、水素4つ・炭素2つ・酸素1つの重さを加えています。

#重合度からセルサイズを決定--------------------------------
print('重合度をどうしますか?(nは4以上)')
#高分子鎖の重合度
n=int(input())
c_leng=(n-1)*3.30
cut_leng=12

#セルサイズ
print('この場合の高分子の最大長さは'+str(c_leng)+'です。セルの一辺の長さはどうしますか?')
print('※今回のカットオフの長さは'+str(cut_leng)+'です。')
cellsize=int(input())
x=cellsize
y=cellsize
z=cellsize

 この部分では粗視化粒子の重合度nとセルのサイズを定義しています。

『n=int(input())』と記述している部分では、プログラムを実行した時にキーボードで入力した数をn代入する、という操作が行われています。こうすることでプログラムを実行したあとで設定を変更できるので便利です!(もし変更する機会が少なければ数を代入すればOKです)

『c_leng』は重合度nに対する高分子鎖の全長を計算しています。(bondポテンシャルの安定距離が3.3Åなので3.3にしています)『cut_leng』では長距離間力の限界距離(カットオフ距離)を代入しています。

後半のセルサイズの部分でも重合度と同じようにその場で入力できるようにしています。ここで注意して欲しいのは『セルサイズとカットオフ距離の関係』です。もしセルサイズを高分子の長さと同じにしてしまうと、周期境界条件にした時に時刻ゼロの時の粒子の距離が近すぎるため過剰な力が加わり、シミュレーションが破綻してしまいます。(これはよくあるミス)そのためセルサイズを決める時は高分子が取りうる最大長さを考慮してセルサイズを決めましょう。

 

 

#各種設定-------------------------------------------------------
#高分子の色んな数
c_bond_n=n-1
c_angle_n=n-2
c_dihed_n=n-3

atoms=n
atom_types=1
bonds=c_bond_n
bond_types=1
angles=c_angle_n
angle_types=1
dihedrals=c_dihed_n
dihedral_types=1

 ここでは粗視化粒子のナンバリングのための変数を記述しています。

最初の3つはボンドポテンシャルや角度ポテンシャルなどの総数を定義しています。

以下の図のように粗視化粒子がn個の時にはボンドポテンシャルがn-1個、アングルがn-2個、となることがわかります。

後半の部分では粗粒子の総数(atoms)、粗粒子の種類(atom_types)というようにlammpsで最初に書かなければならないパラメータの数を定義しています。自分のやりたいシミュレーションに応じて変更する必要があるので、自分で考えてみてください。

print('#各種結合の設定')
print(str(atoms)+' atoms')
print(str(atom_types)+' atom types')
print(str(bonds)+' bonds')
print(str(bond_types)+' bond types')
print(str(angles)+' angles')
・・・

 個々の部分では先ほどまでに設定した変数をターミナルで表示するように命令しています。記述する順番を間違えると、lammpsがエラーを吐くことが多いので注意しましょう。

#変数
#Noの為の変数
no=1
#炭素鎖の奥行方向の座標
c_x0=x/2
#炭素鎖の軸方向の座標
c_y0=y/10
#炭素鎖の高さ方向の座標
c_z0=z*3/4
#粒子のの配置間隔
#wantleng=
#step=n*3.3/wantleng
step=1.0

 ここの部分では高分子の座標を決定するために必要な変数を定義しています。『c_x0』という変数は高分子をセルのどこに配置するかを決める変数です。今回の場合では、『x軸方向はセルの2分の1に配置』『y軸方向はセルの10分の1に配置』『z軸方向はセルの4分の3に配置』というようになっています。セルの境界付近にすると計算が不安定になるので気をつけましょう。

『step』の部分は粗粒子を配置する間隔(Å)を示しています。長すぎても近すぎてもポテンシャルが発散するので気をつけましょう。

#高分子鎖の原子位置
for i in range(n):
    #左の炭素
    print(str(no)+' 1 1'+' '+str(c_x0)+' '+str(c_y0)+' '+str(c_z0)+' 0 0 0')
    no+=1
    c_y0+=step

 ここの部分で高分子の座標を出力しています。

『for i in range(n)』の部分は『変数n(重合度)の回数だけ以下の命令を繰り返す』という意味です。今回は粒子の数が重合度と等しいのでこのような書き方になっています。

このループでは以下のような命令が繰り返されています。

  1. 粒子の位置を記述する
  2. 粒子のナンバリング変数(no)に1を加える
  3. 粒子の位置を変更する

『+=』という記号は右辺に書かれた数を代入するという記号です。原子のナンバリングは1刻みなので1ずつ加算しますが、y座標の位置は任意で決められるので先ほど設定した粗粒子を配置する間隔を代入しています。

 

#Bondの変数
b_no=0
#高分子鎖のBondをつなぐための変数
b_n1=1
b_n2=2


#繰り返し部分
for i in range(c_bond_n):
    b_no+=1
    print(str(b_no)+' 1 '+str(b_n1)+' '+str(b_n2))
    b_n1+=1
    b_n2+=1
print(' \n')

 ここの部分ではボンドポテンシャルを記述しています。

最初の部分ではボンドポテンシャルに必要な数を定義しています。今回の高分子では隣り合う分子にのみbondが発生しているので隣り合う変数を用意します。(b_n1とb_n2)

これらの変数を先ほどと同様にfor i in rangeで繰り返し出力することで高分子の結合が記述できます。

 

ここ以下のコードはボンドでやっていることと変わらないので省略します。

 

さいごに

 いかがでしたか?

コードの解説をしたことがなかったのでざっくりとした説明になってしまいましたが、もっと詳しく教えて欲しいという部分があれば、コメント欄で聞いてください('◇')ゞ

【お勉強】分子動力学法-2(分子動力学法で分かること)

はじめに

今回の記事では、分子動力学法からわかる情報について解説したいと思います!

(次から実際のシミュレーションを再開するので、もう少しお付き合いください!)

 

 

分子動力学法で分かること

直接わかること

まず何も計算しなくても分かることは『原子の位置や速度(一次情報)』です。また速度別に分布をグラフ化すると、平衡状態の速度分布(マクスウェル・ボルツマン分布)も分かります。

 

これらの一次情報から、

静的な物理量:温度, 圧力, 比熱, 動径分布関数など

動的な物理量:拡散係数, 粘性係数など

が分かります!

以下ではこれらの物理量がどのように算出されるのかなどを説明します。

 

・温度

 分子動力学法では『エネルギー等分配則』が成り立ちます。エネルギー等分配則とは、

系の持つ自由度ごとに一定量のエネルギーが配分されるという統計力学の法則(Wikipedia)

 というものです。この法則が成り立つのは、運動エネルギーの平均値を近似し量子力学的な効果が無視できる場合だそうですが、ここでは詳しい話は省きます。

(参考にするならここ?⇒https://eman-physics.net/statistic/equipartition.html)

ここでエネルギー等分配則が成り立つとき、1自由度あたりエネルギーは『1/2kBT』づつ分配されます。つまりn個の粒子がある系では、以下の式が成り立ちます。

f:id:toal_peg_simu:20200114223903p:plain(x,y,z方向に自由度があるため3倍されている)

よって系の温度は粒子の速度から求めることができるんです!
 

・圧力

次に圧力についてですが、圧力は『ビリアル定理』より算出します。ビリアル定理とは、

多粒子系において、粒子が動き得る範囲が有限である場合に、古典力学量子力学系のいずれにおいても成立する以下の関係式のこと

f:id:toal_peg_simu:20200114224607p:plain(Wiki)

 というものです。ここでも複雑な説明は省きますが、運動エネルギーとポテンシャルエネルギーに関係性を与える定理、と思ってもらえればいいと思います。この式を変形すると、

f:id:toal_peg_simu:20200114225309p:plainとなり、f:id:toal_peg_simu:20200114225323p:plainf:id:toal_peg_simu:20200114225331p:plainを代入すると、

f:id:toal_peg_simu:20200114225452p:plain

が算出されます。つまり系の圧力は原子の速度・位置・力から求めることができます。圧力は式を見ればよく分かりますが、ポテンシャル関数に依存しているため圧力の揺らぎは大きくなります。シミュレーションで議論可能な圧力データを得るには長時間の解析が必要になるので大変ですね・・・

 

・動径分布関数

そもそも動径分布関数とは、

等方的な系(または角度依存性を近似的に無視できる系、球対称な系)の中で、ある物理量の分布が原点からの距離 r のみの関数である場合に、その分布を表す関数である。(Wikipedia)

という物理量です。シミュレーションでよく用いられるのは粒子同士の距離の分布です。例えば現実の水分子と水分子のモデルの差異を検証する際に、水素原子と酸素原子の分布で比較したりします。

系全体の動径分布関数g(r)は、

f:id:toal_peg_simu:20200114231011p:plain

で表せます。つまり特定の時刻における粒子の位置から算出することができます。

 

・拡散係数

拡散係数とは、濃度勾配がある場合に単位時間あたり物質が移動する量(流束)を示す物理量です。この物理量が分かると、ある溶媒中の溶質の運動について議論できるようになります。

もし1種類の物質のみから成る系の拡散係数(自己拡散係数)について考えるなら、

f:id:toal_peg_simu:20200114232202p:plain

というアインシュタインの関係式が成り立ちます。つまり各時間の粒子の位置から拡散係数も算出できるんです!

 

周期境界条件

最後に分子動力学法をマクロスケールまで適用するための手法、周期境界条件について説明します。

周期境界条件とは、マクロな系(10^{23}オーダー)をシミュレーションで扱えるミクロな系(10^{9}オーダー)で表せるようにする『賢い』境界の設定です。

f:id:toal_peg_simu:20200114234108p:plain

図で簡単に解説すると上図のようになります。中央のセルのみが計算対象の系で、そのセルの周りに模写した仮想セルを配置します。このセルから粒子が出た場合にはセルの反対側から粒子が侵入するようになります。またセルの境界付近でポテンシャルについて考える時は、仮想セルの粒子も考慮に入れて計算が行われるようになります。

このような条件にすることで有限個の原子からバルクの系を再現できるようになります。ただしポテンシャルのカットオフ距離や系の大きさを適切に定義しなければならなないため、慎重に解析パラメータを設定する必要があるがあります。

さいごに

 いかがでしたか??

今回は分子動力学法から得られる情報とバルクの物性を表現するためのシミュレーションの工夫について解説しました。

疑問点があればコメントお願いします(=゚ω゚)ノ

【お勉強】分子動力学法-1(原子に作用する力)

はじめに(2020/01/12加筆)

今回の記事はlammpsで主に使われる分子動力学法の解説をしたいと思います。

(今回の記事は学問的なものなので問題点があればご指摘頂けると幸いです。)

 

 

 

簡単な説明

分子動力学法(MD,Molecular Dynamics)とは全ての原子の運動について運動方程式を解くことにより、その軌跡を追跡する計算方法です。運動方程式とは『ニュートンの第二法則』の運動方程式です。

f:id:toal_peg_simu:20200106095111p:plain ニュートンの運動方程式 - Wikipedia

 

 この方程式にシミュレーション開始時の『初期位置』と『初期速度』、各時間における『原子に作用する力』を代入すると、次のstepにおける原子の位置と速度を求めることが出来ます。

 

原子に作用する力 is 何?

先ほどの運動方程式に影響する『原子に作用する力』には以下のようなものがあります。

シュレーディンガー方程式(⇒第一原理分子動力学)

f:id:toal_peg_simu:20200106101524p:plain  シュレーディンガー方程式 - Wikipedia

シュレディンガー方程式を解くと電子の波動関数が分かり、そこから『フェルマン・ファイマン力』を求めることが出来ます。この力は原子核に作用するため、原子の運動に影響を与えます。

電子状態が変化するような系(レーザーによる励起など)を取り扱う際に有効です。ただし時間に依存するシュレディンガー方程式は解くことが難しいので、いくつかの拘束条件が必要になります。この拘束条件を系に合わせて設定しなければシミュレーション結果に悪影響が生じます。

 

②ポテンシャル関数(⇒古典分子動力学)

原子に作用する力を『ポテンシャル関数』で近似してシミュレーションをします。有名なポテンシャル関数には『Lennard-Jonesポテンシャル』があります。

U(r)=4\epsilon \left[\left({\frac  {\sigma }{r}}\right)^{{p}}-\left({\frac  {\sigma }{r}}\right)^{{q}}\right] レナード-ジョーンズ・ポテンシャル - Wikipedia

 うまくポテンシャルを設定すればシュレディンガー方程式よりも計算量を減らすことが出来るため、大きな系についても扱う事ができるようになります。一方でポテンシャルを事例に応じて調べる必要があります。

前回の記事でシミュレーションでは古典分子動力学を使用しました。

 

原子の位置と速度を決定する(運動方程式の数値積分)

原子の位置と速度を決定するためには運動方程式に当てはめて解く必要があります。パソコン上で微分方程式を解く場合には差分方程式に変換します。

この差分方程式を解くにはいくつかの手法があります。

Leap frog法

{\displaystyle {\begin{aligned}x_{i+1}&=x_{i}+v_{i}\Delta t+{\frac {1}{2}}a_{i}\Delta t^{2},\\v_{i+1}&=v_{i}+{\frac {1}{2}}(a_{i}+a_{i+1})\Delta t.\end{aligned}}} リープ・フロッグ法 - Wikipedia

二次近似で微分方程式を解く手法である。この手法のメリットは時間可逆性があるということです。また運動のエネルギーが保存されるという性質ももちます。

ただしこの手法では加速度が速度に非依存であり、時間ステップを変えることができないので注意が必要です。

 他にはverlet法などがあります。(解説したかったですが、イマイチ違いが分からないのでまた今度にします)

 

Leap frog法(詳しく)

 Leap frog法の所で『時間可逆性』『エネルギー保存則』のところがイマイチ分からなかった、というコメントをもらったので少し調べてみました。

 

時間可逆性

まず時間可逆性が守られる意味についてなんですが、そもそも現実世界は非・時間可逆性らしいです。分かりやすい資料があったので以下に載せておきます。

https://www.gakushuin.ac.jp/~881791/materials/Irreversiblity09.pdf

https://en.wikipedia.org/wiki/Time_reversibility

ここで注目すべきことは『ニュートンの力学法則で成立している系では、時間の向きが可逆である』ということです。つまり分子動力学法ではニュートンの第二法則が系を支配しているため、時間の可逆性が成り立つことになります。

ここでもう一度数値積分について考えてみましょう。シミュレーションにルンゲクッタ法を用いると、そもそもルンゲクッタ法に時間可逆性がないため、シミュレーション自体も(ニュートン力学の系なのに)時間可逆性が失われてしまいます。このような論理の破綻をなくすために、分子動力学法では時間可逆性がある『Leap frog法』などが使われているのかもしれませんね。

 

エネルギーの保存

修正前の記事では時間可逆性とエネルギーの保存に因果関係があるかのように書いていましたが、実際は独立の性質でした・・・(すみません、間違っていたので書き直しました)

恐らくここで言いたかったのは、各粒子に初期値として与えられた運動エネルギーが長時動シミュレーションをしても変動しないことがleap frog法のメリット、ということが言いたいのではないでしょうか。

 

 

 

さいごに

いかがでしたか。

今回はlammpsでよく使う分子動力学法の簡単な説明と運動の推測手法について説明を行いました。

難しい表現で分からない部分があればコメントなどで指摘してください。

 

 

【lammps】入力ファイルの解説!!_2

はじめに

前回紹介できなかったCONFIGファイルとParamファイルの解説を行います!!

 

 

解説

②CONFIGファイル

#各種結合の設定
100 atoms
1 atom types
99 bonds
1 bond types
98 angles
1 angle types
97 dihedrals
1 dihedral types

 ここではシミュレーションで使用する粒子の数や結合数などを明示しています。

ここで数を間違えるとエラーが起きます。めんどくさいと最初は思っていましたが、原子の位置をプログラミングとかで出力するときに間違いがあるとここで分かります。ありがたいです(笑)

#cell initial size
0 150 xlo xhi
0 150 ylo yhi
0 150 zlo zhi

ここではセルのサイズを定義しています。ここでエラーを吐くことはあまりないですが、セルに粒子を詰めすぎると周期境界条件の時なんかはゼロ距離になって発散するので注意!!

Masses

1 44.0

 ここでは原子の重さを定義しています。最初についている『1』はインデックスです。

Atoms

1 1 1 50.0 0.0 50 0 0 0
2 1 1 50.0 1.0 50 0 0 0
3 1 1 50.0 2.0 50 0 0 0
4 1 1 50.0 3.0 50 0 0 0
5 1 1 50.0 4.0 50 0 0 0
~

 ここでは粒子の位置と粒子の種類を明示しています。

一列目が各粒子のインデックスです。物理量の指定やbondやangleの指定の時に使います。

二列目が分子のインデックスです。今回は高分子が1本しかないので1しかないです。

三列目は粒子の種類です。Massesのところで指定した番号の重さが適用されます。

四列目、五列目、六列目が粒子のxyz軸上の位置を示しています。

それ以降は今回は使わないので無視します。(今回の説明はマニュアルの『compute property』に書かれています)

Bonds

1 1 1 2
2 1 2 3
3 1 3 4
4 1 4 5
5 1 5 6

 ここでは結合している粒子を定義しています。

一列目はBondのインデックスです。

二列目はBondの種類です。Paramファイルで指定したポテンシャルが適用されます。

三列目と四列目が繋がっている粒子のインデックスです。

AnglesとDihedralsの部分も同じように書かれているので省略します!

③Paramファイル

#### References ################################################
pair_style lj/gromacs 9.0 12.0
bond_style harmonic
angle_style cosine/squared
dihedral_style fourier

 ここではシミュレーションで利用するポテンシャル関数を指定しています。自分が使いたいポテンシャルをマニュアルから探して適用します。(各関数の説明は理論の記事で取り上げたいです)

#### PAIR ######################################################
pair_coeff   1 1 0.8066 4.3
#### BOND ######################################################
bond_coeff 1 20.32 3.30
#### ANGLE ######################################################
angle_coeff 1 10.16 130
#### DIHEDRAL ######################################################
dihedral_coeff 1 4 0.4685 1 180 0.04302 2 0 0.07887 3 0 0.02868 4 0

 ここではポテンシャルの関数で使用する係数を指定しています。

ここの係数は論文の書き方やマニュアルによって違うので気を付けましょう。

 

さいごに

いかがでしたか?

ファイルの解説はここで終了です!次回以降の記事は、シミュレーションの基礎知識について書いていきたいと思います(=゚ω゚)ノ

【lammps】入力ファイルの解説!!_1

はじめに

今回は前回の記事で使用したファイルについて解説したいと思います!(※内容が多すぎるので分けました!)

今回解説する時にしたのはlammpdsの公式マニュアル(英語)です。英語読める人はそっちを参考にするのが良いと思うので以下をご覧ください。2000ページあって重いので注意!

https://lammps.sandia.gov/doc/Manual.pdf


 

 

解説

①INファイル

##/bin/sh

 これはシミュレーションに直接関係ないですが、Unix系で使われるシェル「bash」のシンボリックリンク?らしいです。詳しく知りたい人は調べてみてください。

units real
atom_style  molecular

 一行目はシミュレーションで使用する物理量の単位を規定しています。例えばrealスタイルではmassが『grams/mole』となっていますが、siスタイルにすると『kilograms』になります。(単位を間違えると全てがパーになるから気を付けような←)

二行目はシミュレーションで使用する粒子が持つことができる物理量を規定しています。例えばmolecularスタイルでは『only the default values(速度とか原子の種類のコマンド)』を使う事ができます。一方でfullスタイルでは『molecular style + charge』となっているのでmolecularスタイルに加えて、電荷を付与することができます。どれを選べばいいかはマニュアルに書いてあるので『atom_style command』の所を読みましょう。

boundary p p p
read_data   CONFIG
include Param

 一行目ではシミュレーションのセルの境界条件を規定しています。pは周期境界条件を示しています。一文字目がx軸、その次がy,zに対応しています。周期境界条件にすると大きな系を小さなセルで表現することが可能になるので、高分子シミュレーションの場合ではよく使われます。

二行目,三行目は読み込むファイルを示しています。このようにファイルを分けたのは、一つのファイルに全てを記述すると見通しが悪くなるからです。(特に系が大きくなると大変)

neighbor  10.0 bin
neigh_modify  delay 0  every 1  check yes
neigh_modify  page 2000000  one 50000 

ここではポテンシャルを計算するときのパラメータを決定している。 (イマイチよく分かってない部分だが)雑に設定し過ぎると、原子の位置をlostして解析がストップするので注意が必要である。

一行目は『neighbor pageに(cut offに10Åを加えた距離にある)原子までを記録する』ということを示している。

二行目は『各ステップごとに原子の位置を更新する(ただし原子が大きく移動した時のみ)』ということを示している。

三行目は『一つのneighbor pageにpairを2000000個保存して、一つの粒子に対してneighborを最大50000個まで集計する』ということを示している。

(ここでneighbor pageとはポテンシャルを計算をする原子のペアを収納しているファイルではないかと考えられる。)

variable  dt    equal   10
variable  Nstep equal   40000
variable  dNcfg equal   ${Nstep}/100
variable  dNout equal   ${Nstep}/100
variable  Tini  equal   323
variable  Tend  equal   323
variable  Tmass equal   ${dt}*100
variable  Pini  equal   1.0
variable  Pend  equal   1.0
variable  Pmass equal   ${dt}*1000

 ここでは後の計算をする時の変数を設定しています。

dtはシミュレーションの1stepあたりに進める時間、Nstepはシミュレーションをする総step、dNcfgはcfgeファイルを出力するstep間隔、dNoutは物理量を出力するstep間隔、TiniとTendはシミュレーションをする時の温度、Tmassは温度に関するdampingパラメータ、PiniとPendとPmassは圧力の以下略である。

thermo_style  multi
thermo  ${dNout}

 ここでは出力する熱力学量を規定しています。multiにすると、log.lammpsに色んな物理量が記述されます。(詳しくはマニュアル読んで)

二行目では物理量を出力するファイル数を規定しています。

timestep  ${dt}
velocity all create ${Tini} 1 rot yes dist gaussian
fix 1 all nvt  temp ${Tini} ${Tend} ${Tmass}

 ここではシミュレーションをする際の条件を示しています。

一行目は1ステップあたりに進める時間を示しています。大きいほど時間経過が観察しやすいですが、安定しなくなります。(理論の記事でまた説明します)

二行目はシミュレーション開始時に粒子に初速を与えるかどうかを示しています。今回の場合は全ての粒子(all)にTiniケルビンの初速をガウシアン分布で与えています。

三行目はnvt条件(等温等積粒子数変化なし)でシミュレーションを行うことを示しています。

fix engout all print ${dNout} "$(step) $(ke) $(pe)"                   screen no file energy.d
fix strout all print ${dNout} "$(step) $(pxx) $(pyy) $(pzz)"          screen no file stress.d
fix cllout all print ${dNout} "$(step) $(vol) $(lx) $(ly) $(lz)"   screen no file cell.d
dump 3 all cfg    ${dNcfg}   *.cfge mass type xs ys zs mol fx fy fz
#dump 4 all custom ${dNcfg}   *.atom id mol type id xu yu zu fx fy fz
dump_modify 3 element C

ここでは出力するファイルを指定しています。指定した物理量がファイルに出力されます。例えばstress.dの左から3列目はセルのy方向にかかる圧力を示しています。

run    ${Nstep}
write_data CONFIG.OUT nocoeff

ここではシミュレーションの開始と原子位置ファイル(CONFIG.OUT)の出力を指示しています。

 

まとめ(INファイル)

いかがでしたか?

INファイルがシミュレーションの設定をしていることがおわかりいただけたでしょうか?次回はCONFIGとParamファイルの解説をします!!

【lammps】シミュレーション結果を可視化してみよう!

はじめに

今回は前回の記事で作成した出力ファイルを可視化する方法を考えてみましょう!

 

ソフトの紹介

原子構造を可視化する方法は2つあります!

①atomeye

lammpsで排出されたファイルをそのまま可視化することができるソフトでぬるんぬるんに動きます。

(※今回はうまく環境設定ができなかったので飛ばします・・・また今度)

 

②OVITO

同様にlammpsで排出されたファイルをそのまま可視化することができるソフトで、UIがしっかりとしています。

以下のページからインストールすることができます。多くのOS向けに作られているので自分の環境に合ったバージョンを選びます。

https://www.ovito.org/

 

使い方(OVITO)

①OVITOを起動する

こんな感じ。ここんところターミナルしか扱ってないのでUIの優しさが身に沁みますね♡

f:id:toal_peg_simu:20200102213026p:plain

②.cfgeファイルをウインドウにドラッグする

こんな感じです。まじで便利。UI最高。

これは一番最初の粒子の配置です。こんな感じだったんですね。

f:id:toal_peg_simu:20200102213149p:plain

③粒子の色を指定する

粒子の色や粒径を設定しましょう。右の場所でInputのParticle typesでをクリックすると色と粒径が変えられます。

f:id:toal_peg_simu:20200102213552p:plain

 

④画像の保存先を指定する

鉛筆マークの右のボタンを押すとレタリングの設定になります。single frameにするとある瞬間のスナップショットになります。animationにすると連番で写真が作成されます。

f:id:toal_peg_simu:20200102214002p:plain

⑤画像を排出する

 save fileでchooseボタンを押して画像の保存先を指定できる。指定したらRender Active Viewpointを押すと排出が始まります。

ファイルの形式を変えれば動画も出せます!!

 

さいごに

いかがでしたか?前回のシミュレーションでどんな感じになっていたのか分かってもらえたでしょうか(笑)

次回は解析ファイルの解説をしていきたいと思います!!