MD tutorial 3: Ar atoms
I found that the MD simulation using DL-POLY is not the handsoff calculation like a gaussian G03 program.You should makesomething before and after the simulation. But I like it because it
is not blackbox simulation.

2013年 7月 2日 火曜日

1

Ar atom:atomic number 18
electron 1s2 2s2 2p6 3s2 3p6☞ neutral inert gas
atomic weight 39.948

If there is no interaction between Ar atoms,
they are always in the gas state?
Butfrom experimental data

Melting point: 83.8 K, Boiling point 87.3 K
Critical temperature 150.8 K
Enthalpy of vaporization 6.5 kJ mol-1
2013年 7月 2日 火曜日

2

What’s the origin of the interatomic force?

• attractive force: van der Waals force from instantaneous
dipole-dipole interaction

•

•

c
∝− 6
r

in the following you can see the origin of vdW interaction and c is given

3 α1 α2 I1 I2c=
2 (4π�0 )2 I1 + I2

• repulsive force : overlapping of wavefunctions
One can assume that

2013年 7月 2日 火曜日

c�
∝ 12
r
3

2013年 7月 2日 火曜日

4

2013年 7月 2日 火曜日

52013年 7月 2日 火曜日

6

2013年 7月 2日 火曜日

7

2013年 7月 2日 火曜日

8

vdW interaction (or LennardJones) is the one of the mostimportant intermolecularinteraction and itʼs always
used in MD calculation.

2013年 7月 2日 火曜日

9

Lennard-JonesPotential�cc
V (r) = − 6 + 12 = 4�
r
r

�� �
σ 12
r

−

� σ �6 �r

V (σ) = 0
V (2

1/6

Force

F (r)

r0 = 21/6 σ� � ��� σ �5dV (r)σ 11= −= 4� 12(−σr−2 ) − 6(−σr−2 )
dr
r
r
� � �
�
24�
σ 12 � σ �6
= −
2
−
r
r
r

F (r0 ) = −2013年 7月 2日 火曜日

1 1σ) = 4�( − ) = −�,
4 224�
(2/4 − 1/2) = 0,r

F (σ) = −24�/σ
10If there are two Ar atoms
in the system, the equation of
motion becomesOF�iV

r1 - r2

r1

r2

∂V
∂V
∂V
[means (Fi )x = −, (Fi )y = −, (Fi )z = −], i = 1, 2
∂xi
∂yi
∂zi
�
= V (|r1 − r2 |), |r1 − r2 | = (x1 − x2 )2 + (y1 − y2 )2 + (z1 − z2 )2= −∇i V,for example

�∂V∂V∂|r1 − r2 |dV ��1x1 − x2
�
=
=
2(x
−
x
)
=
V
(|r
−
r
|)
1
2
1
2
∂x1
∂|r1 − r2 |∂x1
dr �r=|r1 −r2 | 2|r1 − r2 ||r1 − r2 |Then, equation of motion becomed2
m1 2 r1
dt==

d2
m2 2 r2dt
2013年 7月 2日 火曜日=

x1 − x2
y1 − y2z1 − z2i − V � (|r1 − r2 |)j − V � (|r1 − r2 |)k|r1 − r2 ||r1 − r2 ||r1 − r2 |r1 − r2−V � (|r1 − r2 |)r̂12 , r̂12 =|r1 − r2 |−V � (|r1 − r2 |)V � (|r1 − r2 |)r̂12
11

Lennard-Jones n-l Potential
φ(r) = 4�
φ� (r) =φ(σ) =φ� (r0 ) =φ(r0 ) =d2
mk 2 �rk
dt

2013年 7月 2日 火曜日

�� �
σ n

−

� σ �l �

r
r
� � �
�
�
�
n
l
dφ(r)
�
σ
σ
= −4 n
−l
dr
r
r
r
0
1
� n � n−l0, r0 = σl�� � n
l �
�
�
l n−l
l n−l
4�
−
n
n
�

xk − xiyk − yizk − zi��[φ (rki )i + φ (rki )j + φ (rki )
k]
= −
rki
rki
rki
i
�

12

∂µ

i

ki

V,T

ki

zk − zAr
i atoms:
In the +φcase" of many
(r )
k]
XIV.

PAIR POTENTIALkir

If the interaction between i and j is pairwise, the potential energyki
is given by!
1!
φ(|ri − rj |)φ(|ri − rj |) =
2 i,j
i<j

V=

(103)

The force acting on the atom k can be written as

XV.

LENNARD-JONESPOTENTIAL#$
14%fk = −∇k V = − i ∂ + j ∂ + k ∂
∂xk
∂yk
∂zk
2
2
2 1/2
= [(xi r−+i −(yri j−| yj ) + (zi − zj ) ]j ) |r
ij x=

(104)

V

(105)

∂|r −byr | the1
be described
Lennard-Jones type, for example,= |r − r | 2(x − x )(δ − δ )
(106)
∂x
2
!
"x − x $+ jφ (r ) y#− y $ %#f = − 1[iφ (r ) #
2
rσ 12
r σ 6
φ(r) = 4%z − z−+kφ (r )](δr − δ )
(107)
r
r
!
%
y −y " #x −x#$$i + φ (r )j= −[φ (r )
σ 6
dφ(r)
σ 12
r
r%
"
−
φ (r) = z − z = −24 2+φ (r ) drk]rrr(108)r2!
dXV.
yk − yixk − xi""
LENNARD-JONES POTENTIAL
mk 2 #rk = −i + φ (rki )j[φ (rki )dtby the Lennard-Jones type,
r
r
ki
ki
action can be describedfor example,
i
i

j

k

i

j

−1
"

k

i

j

i

ij

i

ij
"

j

ij

k

ki

φ(r) = 4%

ki

k

"

ki

kj

i

"

ki

i

"

j

kj

ij

i,j
"

ki

ki

i

ij

j

ij

k

i

ki

i

ki

"# $
σ 12

# σ $6 %z −z+φ (r )
k]
r

−
k
i
"r
r
%
" #ki
dφ(r)
%
σ $12 # σ $6"−
φ (r) =
= −24 2
dr
r
r
r ki
2
!
yk − yixk − xi
2013年 7月 2日 火曜日 d#r = −
mi + φ" (r )j[φ" (r )

(109)
(110)

13

c

--periodic boundary condition
axh=ax/2.0d0
ayh=ay/2.0d0
azh=az/2.0d0
do i=1, natom
if (x(i) > axh) x(i)=x(i)-ax
if (x(i) < -axh) x(i)=x(i)+ax
if (y(i) > ayh) y(i)=y(i)-ay
if (y(i) < -ayh) y(i)=y(i)+ay
if (z(i) > azh) z(i)=z(i)-az
if (z(i) < -azh) z(i)=z(i)+az
enddo

Allen and Tildesley, Computer simulation of liquids, Oxford

To neglect the effect of wall, PBC is used.
2013年 7月 2日 火曜日

14When the PBC is set up, the potential should be cut-off,because a molecule in the center cell can interact with the self-image in theneighboring image cell.

To avoid that the total energy is not conserved by introduction of thecutoff, the potential is shifted�v(r) − v(rc ) r ≤ rc
Sv (r) =
0r > rcTo avoid the discontinuity of the force at the cutoff, the force-shifted potential
is sometimes introduced
�
dv(r)
v(r)
−
v(r
)
−
c
dr |rc (r − rc ) r ≤ rc
v SF (r) =0r > rc

The correction by the introduction of cutoff should be done.

2013年 7月 2日 火曜日

15

Allen and Tildesley, Computer simulation of liquids, Oxford

Should consider the interaction between 1(central cell) and 4 (cell E)
2013年 7月 2日 火曜日

16minimum image convention in the program (by MY)
For fortran study I strongly recommend 1) “JIS FORTRAN入門(上）森口繁一

東大出版”

and 2) F90 http://ace.phys.h.kyoto-u.ac.jp/~tomita/education/fortran90/sec0.html

by 富田博之.

cc-------- calculate potential ------------------i-j pair include minimum image -do i=1,natom-1
do j=i+1, natom
xij=x(i)-x(j)yij=y(i)-y(j)zij=z(i)-z(j)if (abs(xij) > rcut)! when |x_i - x_j | > r_cut&
xij=xij-sign(1.0,xij)*ax !
! check the image cell x_i - x_j - sign(xij)*ax
if (abs(xij) > rcut) cycle! sign(1.0,x ) = 1.0 when x > 0, = -1.0 when x < 0
if (abs(yij) > rcut)&
yij=yij-sign(1.0,yij)*ay
if (abs(yij) > rcut) cycle
if (abs(zij) > rcut)&
zij=zij-sign(1.0,zij)*az
if (abs(zij) > rcut) cycle! cycle means skip this do-loop
r2ij=xij*xij+yij*yij+zij*zij
if (r2ij > r2cut) cycle
rij=sqrt(r2ij)
.....
enddo
enddo

2013年 7月 2日 火曜日

17寄り道: Programming のおすすめ again
特にSEで就職する方には必携かもね！

For Fortran study I strongly recommend
1) “JIS FORTRAN入門(上）森口繁一

東大出版”(It takes just one day) and in the second step2) F90 http://ace.phys.h.kyoto-u.ac.jp/~tomita/education/
fortran90/sec0.html

by 富田博之. (this is a kind of dictionary.)•••

FORTRANが生き残っているのには理由がある。：特に数値計算分野
������http://www.komoto.org/fortran/��������http://members.at.infoseek.co.jp/kitaurawa/fortran.html•

標準的なC言語でもよい。��http://www9.plala.or.jp/sgwr-t/index.html �•or2013年 7月 2日 火曜日

速習C言語入門—脳に定着する新メソッドで必ず身につく (単行本) �菅原 朋子 (著)18•

http://www.komoto.org/fortran/�� より

Fortranを使おう「いまさらFortran!?」などと言われそうですが、誤解を恐れずに、あえて書くなら、現在普及している高級言語の中でFortranは最も先進的な言語です(ただし数値計算分野に限る)。
簡単にFortranの歴史を書いておきます。
Fortranの生まれは古く、1957年、IBMのJohn Backusにより世界初の高級言語として開発されました。FORmula
TRANslatorの名が示すように、Fortranは当初から数値計算を目的としていました。
その後、FORTRAN66、FORTRAN77規格ができ、FORTRAN77は数値計算用言語として多くの研究者に愛用されてきました。そのため、数値計算関連の書籍やライブラリでは、FORTRAN77で書かれた物が数多くあります。
最近、FORTRAN77の問題点を改善し、さらに数値計算特有の機能を強化した Fortran90及びFortran95ができました。
ちょっと普及が遅れている気がしますけど。

FORTRAN77
Fortranと言うと、未だにFORTRAN77の印象が強いようです。しかし、ここではっきりと書いておきます。FORTRAN77
は既に時代遅れです。もっとひどく言うなら、「古典だ」と言ってしまって差し支えないでしょう。
根本的に、構造化プログラミングをしたくてもできない言語仕様に問題があります。具体的に問題点を指摘すると、以下のような点が挙られると思います。
* 行番号なしでプログラミングできない* ループブロックを記述できない(goto でループを作るしかない)* 構造体がない
いいかげん、FORTRAN77はやめて、Fortran90を使いましょう。

2013年 7月 2日 火曜日

19

Fortran90/Fortran95
Fortran90では、上記のようなFORTRAN77の欠点がほぼ全て解消されています。それに加えて、数値計
算用にベクトル化や並列化に適した構文が追加されています。どんな特徴があるのか、C/C++と比較してみましょう。
配列演算を書きやすい
たぶん、他の言語と比べた時の、Fortranの最も大きな特徴でしょう。特に、多次元配列を扱おうとすると、他の言語を使うのがばからしくなります。
並列性を明示的に書くことができる
FortranのForAll文や、配列演算文を使うと、コンパイラは演算に並列性があることを簡単に解析できます。そのため、コンパイラは自動ベクトル化や自動並列化を行いやすくなります。
可読性が高い（Cは特に読みにくいのでコメントを多く：MY註）"FORmula TRANslator"と言うぐらいですので、元の数式に近い形でプログラムを書けます。また、
配列演算構文を使うことで、何重ものループを書いたりする必要が無くなり、ループだらけで何をやっているのかよくわからない、といったことにもなりません。
豊富な数値計算用ライブラリ
数値計算に使うライブラリ群は、大抵FortranとC向けに作られています。 NetLib などを見ると、いろいろなライブラリが提供されています。
オブジェクト指向は実装されていない
C++では、オブジェクト指向が使えるようになっています。一方、Fortranでは Fortran2000でオブ
ジェクト指向を実装しようという動きがありますが、まだ標準化されていませんし、対応コンパイラもありません。
オブジェクト指向を使えないのは確かにデメリットですが、そもそも数値計算の場合はオブジェクト指向を使うメリットがあまりないので、あまり気になることは無いように思います。
2013年 7月 2日 火曜日

20

試しに、f(x,y,z) = 2 g2(x,y,z) - 3 g(x,y,z) + 4 を FortranとCで書いてみましょう。
Fortranf(:,:,:) = 2.0*g(:,:,:)**2 - 3.0*g(:,:,:) + 4.0

C
for(z=0; z<N; z++){for(y=0; y<N; y++){for(x=0; x<N; x++){f[z][y][x]=2.0*g[z][y][x]*g[z][y][x] - 3.0*g[z][y][x] + 4.0;}}}
Fortranでは、配列演算文をつかって簡単に書けます。このような簡単な式では、Fortranのメリットは感じにくいかもしれませんが、式が複雑になるにしたがって、Fortranの方が書きやすくなります。
また、Fortranの配列演算文では、各配列要素を並列に演算可能であることを示します。上記の例では、コンパイラはN3の並列性があることを簡単に解析できます。
2013年 7月 2日 火曜日

21•�http://members.at.infoseek.co.jp/kitaurawa/fortran.html�より
ＦＯＲＴＲＡＮの歴史
１９５０年代ＦＯＲＴＲＡＮは，米ＩＢＭ社の John. W. Backus 氏が中心となって，１９５４年から１９５７年にかけＩＢ
Ｍ社のＩＢＭ７０４計算機に対して開発された．ＦＯＲＴＲＡＮという名は，FORmula TRANslating system（数式変換システム）の略称から付けられた．世界初の成功した高級言語と言われている．
当時は，高級言語というものに関して懐疑的な見方が強く，そのようなものは非現実的なものと考えられていた．その根拠として，１）コンパイラ自身が貴重な計算機資源を消費すること，２）生成されるコードの効率性がアセンブリ言語によるものほど期待できないこと，などがあげられていた．しかしながら，ＦＯＲＴＲＡＮはこれらの問題を解決し，すばやく効率的な実行コードを生成することに成功した．この成功以降，高級言語に対する懐疑は消え，幾多の高級言語が登場することになる．盟友ＣＯＢＯＬも１９５０年代末には登場し，異端児ＬＩＳＰも同じ頃ＦＯＲＴＲＡＮによって書かれた．
１９６０年代ＦＯＲＴＲＡＮは使いやすい言語でもあったため利用者に好まれ，ＩＢＭ社以外の計算機にも広く普及し，各社
独自のＦＯＲＴＲＡＮが乱立するに至った．その結果，１９６０年代に言語の統一規格の制定が求めらるようになった．結局ＩＢＭ社のＦＯＲＴＲＡＮ ＩＶを元に，１９７２年，世界初の計算機言語規格ＦＯＲＴＲＡＮ６６が制定された．この結果，ＦＯＲＴＲＡＮの普及はいっそう進んだ．
１９７０年代しかしながら，このころ計算機学科方面から，ＦＯＲＴＲＡＮの言語としての欠陥に対する攻撃がはじまる．１
９６８年には，有名な論文『GO TO Statement Considered Harmful（GO TO 文は有害なり）』が現れ，大論争が巻き起こる．そして『structured programming（構造化プログラミング）』が一世を風靡する流行語となる．また，これ以降『GO TO 文』を否定しＰＡＳＣＡＬを礼賛することで通人ぶる風俗が登場する．
これに対し，ＦＯＲＴＲＡＮ陣営は，GO TO のかわりに『COME FROM 文』を提案してみたり，冗談で対抗した．１９８０年に制定されたＦＯＲＴＲＡＮ７７規格は，プログラム制御に関しては BLOCK IF 文すなわち『IF...THEN...ELSE...END IF』という構造化命令だけを取り入れた冗談のような規格である．これだけでは『GOTO 文』を排除できないので，これ以降も『GO TO 文』を満喫する事に成功する．

2013年 7月 2日 火曜日

22

１９８０年代
ところが，そうこうするうち計算機学科ではＰＡＳＣＡＬ対Ｃの内ゲバがはじまっており，ＦＯＲＴＲＡＮ陣営が『RealProgrammers Don't Use PASCAL（本物のプログラマーはＰＡＳＣＡＬを使わない）』などとほざいているうちに，いつの間にかＦＯＲＴＲＡＮは相手にされないようになっていた．１９８１年にはＣ言語のBrian W. Kernighan が『Why Pascal is Not My
Favorite Programming Language（なぜＰＡＳＣＡＬはわたしの好きなプログラミング言語じゃないのか）』を出す．これ以降
Ｃ派がクーデターに成功し，世界がＵＮＩＸ ＷＳ（ワークステーション）とＣまみれになって行くことが運命づけられる．
しかしながら，ＦＯＲＴＲＡＮ陣営は世間の流れにも気づかず，小文字を軽蔑しつつメインフレーム上での楽しい世界にうつつを抜かしていた．さらに１９８０年代のスーパーコンピュータはベクトル計算機であり，この世に存在するベクトル化コンパイラーはＦ
ＯＲＴＲＡＮだけであったので，何の不安もなく田舎地主のごとく気ままに暮らしていられた．
一方，このころＦＯＲＴＲＡＮ７７の中途半端さを憂い，新しい言語規格ＦＯＲＴＲＡＮ８ｘ制定に向けての会議が開かれるが，守旧派対革新派，ＩＢＭ対ＤＥＣ，その他諸々の対立軸でもめにもめ，ずるずると時は流れ，いつの間にか１９８０年代は終わってしまう．結局，１９９１年になって自民党的な素敵な解決策によって新規格ＦＯＲＴＲＡＮ９０がまとめ上げられることになる．
すなわち，革新派にしたがって，構造化プログラミングの必要要素をすべて取り入れると同時に，守旧派にしたがって，これまでのＦＯＲＴＲＡＮ７７規格をすべて包括することになった．
１９９０年代ＦＯＲＴＲＡＮはＦＯＲＴＲＡＮ９０となって Module 概念を取り入れた構造化言語へと進化し１９７０年代の言語水準に到達した．ＦＯＲＴＲＡＮ９０は，理念的だった初期の構造化言語と異なり，様々な実地経験に基づいて規格され，ＦＯＲＴＲＡＮ利用者の様々な不満が解決されるように拡張された．言語の抽象度は上がり，概念をコード化するのが以前より容易になった．しかし，守旧派の多いＦＯＲＴＲＡＮ陣営では未だＦＯＲＴＲＡＮ６６に固執するものいて，遅々として普及は進まなかった．（但しＦＯＲＴＲＡＮ６６は，ほぼＦＯＲＴＲＡＮ７７に包括され，ＦＯＲＴＲＡＮ７７は完全にＦＯＲＴＲＡＮ９０に包括されるので，ＦＯＲＴＲＡＮ６６プログラムは，すなわちほぼＦＯＲＴＲＡＮ９０プログラムであるが．）１９９７年には，ＦＯＲＴＲＡ
Ｎ９０の規格の不備を埋める Bug Fix 的規格としてＦＯＲＴＲＡＮ９５規格が制定され，より完成度を高めている．
２０００年代ＦＯＲＴＲＡＮは，今後もしぶとく生き残ると見られている．その理由は，ＦＯＲＴＲＡＮ使用者の間ではプログラマが即そのプ
ログラムのユーザであることが多い，数学概念との親和性が高い，配列操作に長けている，１９５０年代のコードが未だに残るほど物持ちが良い，ユーザに守旧派思考が強く他の言語やプログラミング概念を知らないし知ろうともしない，などがあげられている．
現在ＦＯＲＴＲＡＮは，オブジェクト指向の概念を取り入れたＦＯＲＴＲＡＮ２０００規格を制定中である．
2013年 7月 2日 火曜日

23

閑話休題

LJ systems
For Ar atoms, m = 6.63x10-26 kg, ε / kB = 119.8 K, andσ = 3.405x10-10 m.

For Ne atoms, m = 3.35x10-26 kg, ε / kB = 47.0 K, andσ = 2.72x10-10 m.

For Kr atoms, m = 13.92x10-26 kg, ε / kB = 164.0 K, andσ = 3.83x10-10 m.

kB = 1.3806503x10-23 J K-1

In the usual LJ systems, one can use the reduced units,
because the system is universal.
In the DL-POLY the following units are used.

2013年 7月 2日 火曜日

24

2013年 7月 2日 火曜日

25

Lennard-Jones system: phase diagram
+

+

+

+

Molecular Dynamics Simulation, J. M. Haile, John-Wiley 1992
2013年 7月 2日 火曜日

26

LJ system is universal:
two parametersρσ 3 , kB T /�ρσ 3 = 0.01, kB T /� = 1.0
ρσ 3 = 1.1, kB T /� = 1.0
ρσ 3 = 0.8, kB T /� = 1.0

2013年 7月 2日 火曜日

27

initial positions of Ar atoms• First we should determine the numbers of atoms
in the cell. (e.g. 864 atoms in liquid and solid
phase.)• In gas phase the density is low, you can put the Ar
atoms randomly using random number
generator.* (* The routine is very important, you should use
good one in Wall-Street if you don’t want to lose your money.)• In solid phase it is known that they make a facecentered cubic(fcc) lattice. So we should put Ar
atoms on the fcc lattice points. It is easy to dothat with your own program.

• In liquid phase we can start from the fcc solid.
2013年 7月 2日 火曜日

28

fcc lattice generation Fortran program cc
by MY (I hope you can understand the language.)
Of course you can make it by hand if you like.(laughter)

c
c

c

c

c

--- fcc lattice generation Fortran program by MY# of atoms:natom
nx=6
ny=6
nz=6
natom=nx*ny*nz*4
write (6,*) 'natom=', natom--- volume
sigma=3.405d0
rho_real= 1.1d0/sigma/sigma/sigma
volume_real=real(natom)/rho_real
pow=1.0d0/3.0d0
box_real=volume_real**pow
ax=box_real
ay=box_real
az=box_real
volume=ax*ay*az
write (6,*) 'ax,ay,az,volume=',ax,ay,az,volume
half of box size
axh=ax/2.0d0
ayh=ay/2.0d0
azh=az/2.0d0
primitive cubic cell
box_nx=box_real/real(nx)box_ny=box_real/real(ny)box_nz=box_real/real(nz)2013年 7月 2日 火曜日

cc---- 4 atoms in cubic cell: fcc(0,0,0), (1/2,1/2,0), (1/2,0,1/2), (0,1/2,1/2,1/2)
do i=1, nx
do j=1, ny
do k=1, nz
icount=icount+1x(icount)=real(i-1)
y(icount)=real(j-1)
z(icount)=real(k-1)
write (6,*) icount,x(icount),y(icount),z(icount)
icount=icount+1x(icount)=real(i-1)+0.5d0y(icount)=real(j-1)+0.5d0z(icount)=real(k-1)
write (6,*) icount,x(icount),y(icount),z(icount)
icount=icount+1x(icount)=real(i-1)+0.5d0y(icount)=real(j-1)
z(icount)=real(k-1)+0.5d0
write (6,*) icount,x(icount),y(icount),z(icount)
icount=icount+1x(icount)=real(i-1)
y(icount)=real(j-1)+0.5d0z(icount)=real(k-1)+0.5d0
write (6,*) icount,x(icount),y(icount),z(icount)
enddo
enddo
enddo
if (icount .ne. natom) stop 999-------- fcc lattice --box shift to -axh to axh and so on.
do i=1, natom
x(i)=x(i)*box_nx-axh
y(i)=y(i)*box_ny-ayhz(i)=z(i)*box_nz-azh
write (6,'(3f15.7)') x(i),y(i),z(i)
enddo
29

• After the atomic

positions are given, you
can calculate the forces
if you know the interatomic potential.

Y

Activefcc lattice that
contains 864 atoms

Z• And if you can integrate
the equation of motion
numerically, you cancalculate the trajectory.

• But we should define

the velocity (If we selectan option, we do not needvelocity and force. But I think
it is better to prepare themin this simplest case.)
2013年 7月 2日 火曜日

30

velocity (momentum)
temperature of the system
N
N
2
�
�
mi vi
p2i
3N kB T
�
�=�
�=
2
2mi
2
i=1
i=1

Maxwell-Boltzmann distribution
�vi = (vix , viy , viz )
�
�1/2
�
�
2
mi
mi vixρ(vix ) =exp −2πkB T
2kB T

viy , viz is the same

momentum of the system → 0
N
�

p�i =

i=1

2013年 7月 2日 火曜日

N
�
i=1

mi�vi = 0c------ MD simulation parameter ---temp=300.0d0! 300 K
bk=1.3806503d-23!kb J K-1
war= 6.63d-26!m_Ar kgc-------- initial velocity ----------c
M.B. Distrubution function
c gasdev is random number generatorc Normal distribution:mean=0,dispersion=1
c 1/(sigma)/(2pi)**0.5 exp(-(x-<x>)**2/2/sigma**2)
cc=sqrt(bk * temp / war)
sumx=0.0d0
sumy=0.0d0
sumz=0.0d0
do i=1, natom
vxx=gasdev(idum)
vyy=gasdev(idum)
vzz=gasdev(idum)vx(i)=vxx*ccvy(i)=vyy*cc
vz(i)=vzz*cc
write (6,*) vx(i),vy(i),vz(i)sumx=sumx+vx(i)sumy=sumy+vy(i)sumz=sumz+vz(i)enddoc---- total momentum => zero----------c
psum is kinetic energy of the system
psum=0.0d0
sumx=sumx/real(natom)
sumy=sumy/real(natom)
sumz=sumz/real(natom)
do i=1, natomvx(i)=vx(i)-sumx
vy(i)=vy(i)-sumyvz(i)=vz(i)-sumz
psum=psum+war*(vx(i)**2+vy(i)**2+vz(i)**2)/2.0d0
enddo
tempinst=2.0d0*psum/3.0d0/real(natom)/bk
write (6,*) psum,tempinstcwrite (6,*) itime
write (6,*) psum,tempinst31

c-------------------------------------------c gasdev is random number generatorc Normal distribution:mean=0,dispersion=1c----------------------------------------------FUNCTION gasdev(idum)
INTEGER idum
REAL*8 gasdev
CU USES ran1
Cran2 use ran2
INTEGER iset
Cran2
REAL*8 fac,gset,rsq,v1,v2,ran1
REAL*8 fac,gset,rsq,v1,v2,ran2
SAVE iset,gset
DATA iset/0/
if (iset.eq.0) then
Cran2 1v1=2.*ran1(idum)-1.
Cran2
v2=2.*ran1(idum)-1.
1
v1=2.*ran2(idum)-1.
v2=2.*ran2(idum)-1.
rsq=v1**2+v2**2
if(rsq.ge.1..or.rsq.eq.0.)goto 1
fac=sqrt(-2.*log(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0endif
returnEND

You can get these codesfrom “Numerical recipe” by Press et al.
http://www.nr.com/

The above codes are modified by MY.
2013年 7月 2日 火曜日

c--------------------------------------------cran2 is uniform randum number generatorcbetween 0 and 1.
c ran2 random deviate by L'Ecuyer long period plus shufflec---------------------------------------------------FUNCTION ran2(idum)
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIVREAL*8 ran2,AM,EPS,RNMXPARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,*IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,*NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) thenidum=max(-idum,1)
idum2=idum
do 11 j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
if (j.le.NTAB) iv(j)=idum
11
continue
iy=iv(1)endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
returnEND
32

random number generator check: ran2
1000 samples

60
50

I(x)

40
30
20
10
0
0.0

0.1

0.2

0.3

0.4

0.5

0.6

x

500

0.7

0.8

0.9

1.0

I(x)

400
300

10000 samples

200
100
0
0.0

0.2

0.4

0.6
x

0.8

1.0

5000
4000

I(x)

100000 samples
3000
2000
1000
0
0.0

0.2

0.4

0.6

0.8

1.0

x

2013年 7月 2日 火曜日

33

ρ(vix ), ρ(viy ), ρ(viz )
40

v̲x
v̲y
v̲z
MB

300 K

6
5

30

-9

3

20

x10

frequency

4comparison between thevelocity distribution
from random number
generator and
M-B distribution

2
10
1
0
-400

-200

0

200

400

ρ(�v )dvx dvy dvz → 4πv 2 ρ(v)dv

v / ms-1

velocity distribution (MD)MB distribution 300 K

50-3

3.0x102.5

40= 10 Å × 10
10

−12

ps

−1

= 0.01 Å ps−1fequency

1ms

−1

30

1.5

20

1.0

300 K

10

MB distribition

2.0

0.5

0

0.0
0

200

400

600

800

1000

v / ms-1

2013年 7月 2日 火曜日

34

Equation of motion: how to solve it
Now we know{�ri }, {�vi }, {F�i } at time t.
How can we know{�ri }, {�vi }, {F�i } at time t + δt.
If the classical trajectory is continuous, the estimate of the positions at t ± δt may be
given by�ri (t + δt)�ri (t − δt)(δt)3 �(δt)2
�ai (t) +bi (t)...
= �ri (t) + δt�vi (t) +2
6(δt)3 �(δt)2
�ai (t) −bi (t)...
= �ri (t) − δt�vi (t) +2
6

By adding and subtracting, we can get�ri (t + δt) + �ri (t − δt)=�ri (t + δt) − �ri (t − δt)=2�ri (t) + (δt)2�ai (t) + O[(δt)4 ](δt)3 �bi (t)
2δt�vi (t) +3Verlet equtation (not recommended)�ri (t + δt)==

2013年 7月 2日 火曜日

2�ri (t) − �ri (t − δt) + (δt)2�ai (t) + O[(δt)4 ]�
2 fi (t)+ O[(δt)4 ]2�ri (t) − �ri (t − δt) + (δt)
mi
35The velocities do not appear in the above formula. The velocities are not needed to computethe trajectories, but they are useful for estimating kinetic energy. They may be obtained by

�vi (t) =�ri (t + δt) − �ri (t − δt)+ O[(δt)2 ]2δtLeapfrog (LF) also used in the electromagnetic wave simulation FDTD Yee algorithm
In this method the mid-step velocity and the positions are updated separately. Using the Taylor
expansionsδt�vi (t − )2δt�vi (t + )2δtδt�vi (t + ) − �vi (t − )
2
2
δt
δt
�vi (t + ) + �vi (t − )22δt(δt)2 �= �vi (t) − �ai (t) +bi (t) + O[(δt)3 ]2
8δt(δt)2 �= �vi (t) + �ai (t) +bi (t) − O[(δt)3 ]2
8= δt�ai (t) + O[(δt)3 ](δt)2 �= 2�vi (t) +bi (t)
4(δt)2(δt)3 ��ai (t) +�ri (t + δt) = �ri (t) + δt�vi (t) +bi (t)...
2
6

2013年 7月 2日 火曜日

36

Then we have�fi (t)δtδtδt3
�vi (t + ) = �vi (t − ) + δt�ai (t) + O[(δt) ] = �vi (t − ) + δt
+ O[(δt)3 ]
2
2
2
mi
�
�
δt
(δt)2 ��ri (t + δt) = �ri (t) + δt �vi (t) + �ai (t) +bi (t) +O[(δt)4 ]
2
6
��
�
�
�
��
vi (t+ δtvi (t)+ δt2 )=�
2 ai (t)+(δt)2 �8 bi (t)δt) + O[(δt)4 ]2δt��vi (t + δt
)
+
v
(t
−
i
2
2)
+ O[(δt)2 ]2= �ri (t) + δt�vi (t +�vi (t) =The calculation procedure is��vi (t − δt/2)
�fi (t)2013年 7月 2日 火曜日

�⇒ �vi (t + δt/2),��ri (t)�vi (t + δ/2)�⇒ �ri (t + δt),�

�vi (t − δ/2)
�vi (t + δ/2)�⇒ �vi (t)37

Velocity Velret!!
The velocity Velret method is recommended because of the memory storage, its numerical stability.
First, the new velocities at mid-step time t + δt/2 are calculated in the following way.
δtδt �fi (t)δt�vi (t + ) = �vi (t) + �ai (t) = �vi (t) ++ O[(δt)2 ]
2
2
2 mi

The new positions are calculated using these velocities at the mid-time-step.
(δt)2�ri (t + δt) = �ri (t) + δt�vi (t) +�ai (t) + O[(δt)3 ]2δt= �ri (t) + δt�vi (t + )2��δt �fi (t)= �ri (t) + δt �vi (t) ++ O[(δt)3 ]2 mi
The forces and accelerations (�f i = mi�ai ) at time t + δt are then computed from �ri (t + δt), andthe velocity at time t + δt is estimated byδtδt�vi (t + δt) = �vi (t + ) + �ai (t + δt)2
2δtδt �fi (t + δt)= �vi (t + ) +
2
2
mi
δt �
= �vi (t) +[fi (t + δt) + �fi (t)] + O[(δt)2 ]2mi
2013年 7月 2日 火曜日

38

The calculation procedure is




�
v
(t)
�
(t)
r
i
 i



�fi (t)�vi (t)⇒ �ri (t + δt) =⇒ �fi (t + δt) →⇒ �vi (t + δt)
 �

 �

fi (t)fi (t + δ)

To avoid the round off error , small number amounts
should be added separately.
�ri (t + nδt) = �ri (t) + δt

n
�

k=1

n

(δt)2 � �
�vi (t + (k − 1)δt) +fi (t + (k − 1)δt)
2mi
k=1

n
�
δt � ��
�
�vi (t + nδt) = �vi (t) +fi (t + kδt) + fi (t + (k − 1)δt)
2mi
k=1

2013年 7月 2日 火曜日

39

Ar atoms MD simulation input files• The followings are the examples.
• I think it is very important to “do it yourself”

microcanonical ensemble NVE

2013年 7月 2日 火曜日

40

1.0

potential / kJ mol-1

0.5

0.0

-0.5

-1.0
5

10
r / Angstrom

15

20

20

2013年 7月 2日 火曜日

10
potential / kJ mol-1

potential min. /100=> r = 10 Angstrom
potential min. /1000=> r = 15 Angstrom=> rcut
rcut = 8 Å is enough

0

-10

-3

-20x10

5

10

15

20

r / Angstrom

41

CONTROL

LJ Argonintegrator velocity verlet
temperature
120.00ensemble nve
steps
5000
equilibration
1000
multiple step
1
scale
10
print
50
stack
100
stats
10
rdf
1
timestep0.004312
cutoff8.00
delr width
0.5000rvdw cutoff 8.00
no electrostatics
shake tolerance 1.0E-5
quaternion tolerance 1.0E-5
print rdf
job timeclose time

600.00100.00finish2013年 7月 2日 火曜日

42

FIELD
LJ Ar
units kJ
molecular types 1
Argon
nummols 864
atoms 1
Ar
39.9480.0
finishvdw 1
ArAr lj 0.99607
close

2013年 7月 2日 火曜日

1

3.405

43

CONFIG2013年 7月 2日 火曜日

Ar LJ system solid 120K rho1.12131.4164760186
0.0000000000
0.0000000000
31.4164760186
0.0000000000
0.0000000000
Ar 1-15.7082380093-15.7082380093-1.4841599823-1.3546851956
0.0000000000
0.0000000000
Ar 2-13.0901983411-13.0901983411-1.1985052697
1.1185526954
0.0000000000
0.0000000000
Ar 3-13.0901983411-15.7082380093-0.2902447339-1.6010119044
0.0000000000
0.0000000000
Ar 4-15.7082380093-13.0901983411
0.2025070266
2.6507153672
0.0000000000
0.0000000000
Ar 5-15.7082380093-15.7082380093
1.6752475218
3.0135226089
0.0000000000
0.0000000000
Ar 6-13.0901983411-13.0901983411-1.8792819470
0.5181167823
0.0000000000
0.0000000000
Ar 7-13.0901983411-15.7082380093-1.9572185128-2.5524163519
0.0000000000
0.0000000000
Ar 8-15.7082380093-13.0901983411-3.0249123171
0.5990315258
0.0000000000
0.0000000000
Ar 9-15.7082380093-15.7082380093-1.1606271365
1.9297868123
0.0000000000
0.0000000000
Ar 10-13.0901983411-13.0901983411
1.2247598317
1.0607088519
0.0000000000
0.0000000000

0.0000000000
0.0000000000
31.4164760186-15.7082380093
0.3078685086
0.0000000000-15.7082380093-0.1253282262
0.0000000000-13.0901983411-1.5683169323
0.0000000000-13.0901983411-1.8006851740
0.0000000000-10.4721586729-2.3565592798
0.0000000000-10.4721586729-2.2284424691
0.0000000000-7.8541190046-1.6063832355
0.0000000000-7.8541190046-1.6698931477
0.0000000000-5.2360793364-0.8672368852
0.0000000000-5.2360793364
0.3416609074
0.0000000000

44

Now we can playaround!!

2013年 7月 2日 火曜日

45

Exercise0) The three input files and input-file-make-fortran source
codes DLPOLY_LJinput.f
can be downloaded fromhttp://www.chem.konan-u.ac.jp/applphys/DL_POLY/

2013年 7月 2日 火曜日

46

From the phase diagram
1) DO the MD simulation of Ar solid ρσ3 = 1.1
at 120 K and get radial distribution function(print rdf option in CONTROL file).
2) DO the MD simulation of Ar liquid ρσ3 = 0.8
at 120 K and get the radial distribution function.
3) DO the MD simulation of Ar gas ρσ3 = 0.01
at 120 K and get the radial distribution function.
4) DO the MD simulation of Ar super criticalstate ρσ3 = 0.3 at 180 K and get the radial
distribution function.
2013年 7月 2日 火曜日

47

Structure
Ar solidρσ3 = 1.1 at 120 K

4

2.0

3

1.5
g(r)

g (r)

Ar liquidρσ3 = 0.8 at 120 K

2.5

2

1.0

1
0.5

0

0.0

0

2

4

6

10

12

14

0

2.0

2

4

6

8

10

12

14

r / Angstrom

Ar gasρσ3 = 0.01 at 120 K
average may not be enough.

2.5

1.5

Ar super critical state
ρσ3 = 0.3 at 180 K

1.5

g(r)

g(r)

8
r / Angstrom

The scs is between
gas and liquid.

1.0

1.0
0.5

0.5
0.0

0.0
0

2

4

6

8

10

12

14

0

r /Angstrom

2

4

6

8

10

12

14

r / Angstrom

The pair distribution function g(r ) gives a probability of finding a pair of atoms a distance r apart, relative to the probability
expected for a completely random distribution at the same�density.
���Vδ(r − rij )g(r) =
N2i
j�=i

2013年 7月 2日 火曜日

48

Dynamics: root of mean square displacement
The mean-square-displacement(MSD) is defined by�� ��|ri (t + t0 ) − ri (t0 )|2 ≡ |ri (t) − ri (0)|2The MSD may be independent on the t0 in the thermal equilibrium condition.
For the short t,ri (t + t0 ) − ri (t0 ) � v(t0 )t��
3kB T 2|ri (t) − ri (0)|2 � �vi2 (0)�t2 =tm

In the longer time limit, the MSD becomes�|ri (t) − ri (0)|

2

�

� 6Dt

In the NVE-MD calculation the center of the mass of the system is fixed, so each particles do not diffuse
freely. By this constraint the MSD deviates from the above linear relation at the longer t limit.
The diffusion constant is also given by

1
D=
3

� ∞
0

dt�vi (t) · vi (0)�In practice these averages would be computed for each of the N particles in the simulation,
the results added together, and divided by N. Note that in the computation of MSD, it isimportant to use the trajectory which have not been subjected to PBC.
2013年 7月 2日 火曜日

49

Be careful that MSD should be
linear with time!! In thesecase no checks were done!!

In the OUTPUT you can find,the diffusion constant is given.

Ar solid ρσ = 1.1 at 120 K
Approximate 3D Diffusion coefficients (10^-9 m^2 / s)
0.0000E+00
3

Ar liquid ρσ = 0.8 at 120 K
Approximate 3D Diffusion coefficients (10^-9 m^2 / s)
4.0951E+00
Ar gas ρσ = 0.01 at 120 K
Approximate 3D Diffusion coefficients (10^-9 m^2 / s)
4.2605E+02
3

3

Ar supercritial state ρσ = 0.3 at 180 K
Approximate 3D Diffusion coefficients (10^-9 m^2 / s)
3.1635E+01
3

2013年 7月 2日 火曜日

50

structure and dynamics
to see the Movie: I like the sofware such as
Jmol http://jmol.sourceforge.net
and/or
VMD http://www.ks.uiuc.edu/Research/vmd/
and If we use VMD, it can read HITORY file
directory.

2013年 7月 2日 火曜日

51

Solid ρσ3 = 1.1 at 120 K

Gas. ρσ3 = 0.01 at 120 K
It is like “Hotaru” firefly2013年 7月 2日 火曜日

Snapshotliquid ρσ = 0.8 at 120 K
3supercritical state ρσ = 0.3 at 180 K
352

Animations (so what's?)
VMD (visual MD) is used here.
HISTORY file can be easily visualised.

Please visit

http://www.ks.uiuc.edu/Research/vmd/
and can export animated GIF by imagemagik!

2013年 7月 2日 火曜日

53

animation2013年 7月 2日 火曜日

Gas. ρσ3 = 0.01 at 120 K
It is like “Hotaru” firefly54

animation2013年 7月 2日 火曜日

Gas. ρσ3 = 0.01 at 120 K
It is like “Hotaru” firefly54liquid ρσ = 0.8 at 120 K
3

2013年 7月 2日 火曜日

55

liquid ρσ = 0.8 at 120 K
3

2013年 7月 2日 火曜日

55

Solid ρσ3 = 1.1 at 120 K

2013年 7月 2日 火曜日

56

Solid ρσ3 = 1.1 at 120 K

2013年 7月 2日 火曜日

56

supercritical state ρσ = 0.3 at 180 K
3

2013年 7月 2日 火曜日

57

supercritical state ρσ = 0.3 at 180 K
3

2013年 7月 2日 火曜日

57

This is the end.
I think it takes time to make a input files if
you make by your-self even in this simple
system.

• In the next tutorial I will try to include longrange coulomb interaction. I will do a hightemperature molten salts (like LiCl and/or
NaCl) MD simulatiuon.

2013年 7月 2日 火曜日

58