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

はじめに

皆さん、こんにちは！

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

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

pythonとは？

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

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

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

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

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

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

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

(※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ずつ加算しますが、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で繰り返し出力することで高分子の結合が記述できます。

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

さいごに

 いかがでしたか？

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