■
ある種の位相的場の理論(Cohomological Field Theory)の"分配関数"(数学では、descendant potentialという呼び方をする時もある)は
(1)何らかの可積分階層のτ関数になっている
(2)Virasoro constraintsや、より一般にW-constraintsと呼ばれる(通常、線形偏微分方程式系で書かれる)条件を満たす
(3)対数は"free energy"で、その無分散極限はWDVV方程式を満たす(prepotential)
(4)"free energy"部分は、種数展開が出来て、それぞれ種数gのGromov-Witten不変量の母関数になっている(種数0部分はprepotential)
という性質を持つということが経験的に分かってるっぽい。幾何学的にはGW不変量の母関数である一方、代数的にみると、Virasoro代数やW代数は、可積分階層のnon-isospectral symmetryとして現れ、Virasoro/W-constraintsは、"大体"最高ウェイト条件となっている。全く異なる由来を持つ2つの量が一致するということで、面白いのだけど、そうなる理由は分からない
具体例としては(予想も含めると)
・Kontsevich-Witten tau function(二次元量子重力理論の分配関数らしい。KdV階層の解)
・Hurwitz tau function(数学的に発見されて、物理的な模型の分配関数として作られたわけではないっぽい。KP階層の解)
・double Hurwitz numberの母関数(戸田階層の解)
・位相的極小模型の分配関数(2D topological minimal matter+topological gravityと解釈できるらしい。Gelfand-Dikii階層/Drinfeld-Sokolov階層)
・位相的シグマ模型の分配関数(よく知らないのだけど、計算できる例は限られていて、射影空間とかtoric Calabi-Yauとかだと分かるっぽい?旗多様体は?)
・様々な行列模型の分配関数
・Nekrasov分配関数?(まぁ、全然知らない)
・2D large N pure SU(N) Yang-Mills理論の分配関数?(Migdalの公式は、無限個の時間変数を含まないので、何か"変形"が必要)
あたり。必ずしも、"位相的"とは思われてなかった理論も含まれる。
1990年代以降、"数理物理"(ミラー対称性とかSW理論)は、幾何学の言葉で書かれることが増えて、1980年代の話題(ソリトン,CFT,可解格子模型)と断絶してしまった感がある(CFTのことを知らない幾何学の人も結構いるらしい)けど、ちゃんと連続した話ではある
コヒーレント状態の有限群での類似物
1970年頃、PerelomovやGilmoreたちが、一般のLie群に対する、generalized coherent stateという概念を考えたのだけど、有限群の場合にも、コヒーレント状態の類似物を考えられるっぽい
参考)coherent state物語
http://d.hatena.ne.jp/m-a-o/20130303#p2
v_1 = (1,0), v_2 = (-1/2 , √3/2) , v_3 = (-1/2 , -√3/2)と、w=(x,y)に対して
w= c_1 v_1 + c_2 v_2 + c_3 v_3
の解(c_1,c_2,c_3)を考える。この解は無数にあるけど、擬似逆行列によって、"自然な"解
がある。この解は
とも書ける。こっちの方が、計算量的にも計算精度的にも優れてるけど、普通は、こういう関係は成立しなくて、こういう風に書ける理由は
という式にある(擬似逆行列の言葉でみれば、転置行列が定数倍を除いて、擬似逆行列に一致するという条件と同じ)。これが、調和振動子のコヒーレント状態とかで、よくみかける式
の類似物(resolution of the identityとか呼ばれる式)
Lie群のコヒーレント状態では、Lie群が、コヒーレント状態の集合に推移的に作用していて、過剰完全基底をなしていた。今の場合、v_1,v_2,v_3は、平面上の120度回転で互いに移りあう。Lie群の場合は、既約ユニタリ表現の最高ウェイトベクトルに群を作用させて、過剰完全基底を得ていたけども、有限群の表現では、最高ウェイトベクトルに相当するものがなく、特に、どのベクトルを選んでもいいらしい(連続ウェーブレット変換のマザーウェーブレットと似たような状況)
有限群の表現は、何でもいいというわけではなくて、既約表現である必要がある。上の3次巡回群の表現は、複素数体上では既約でないけども、実数体上では既約になっている。有限群Gの(適当な体上の)既約表現(ρ,V)で、表現空間には内積(複素表現の場合にはユニタリ内積)が入ってるとして、規格化されたベクトルvを取って、軌道G(v)={ρ(g)v | g ∈ G}を考える。この時、G(v)は、Vの生成系となっている。何故なら、G(v)が生成するベクトル空間Wは、Vの不変部分空間であるけども、(ρ,V)の既約性より、Wは{0}かV自身で、vは0でないので、W=Vとなる。既約でないためにダメな例として、2次巡回群が、xy平面上に鏡映(x,y)->(x,-y)によって作用しているとすると、vの選び方によっては、G(v)が生成する不変部分空間は、一次元になってしまう
G(v)自体は、複数回同じ点を通る場合もあるけど、少し試すと、有限群Gと有限次元既約表現(ρ,V)について、vを長さ1のベクトルとすると、任意のベクトルwに対して
が成り立つという予想が立つ。これは、"有限ウェーブレット変換"とか呼んでもいいような代物だと思う
もう一つ、殆ど同じことではあるけど、wに対して、係数を対応させる上の計算が、一般に擬似逆行列による計算と一致することを言いたい。これは
を示せば、両方一気に言える(e_iは正規直交基底)
証明のアウトラインは以下の通り
が、VからVへの恒等写像の定数倍であるということは、intertwinerであるということなので、逆にこいつがintertwinerであることを言えば、(ρ,V)の既約性とSchurの補題から定数倍の因子を除いて、恒等写像に一致することは分かる。けど、これだと定数倍の因子が出ない。今、vを正規直交基底の一つとすると、この作用素は、vを別の基底ベクトルに取り換えても同じになるので、
となって、 定数倍の因子まで決定できる
Gabor変換や連続Wavelet変換は、最小2-normを与えるという性質で特徴づけられ、擬似逆行列は(解が無数にある不定系に対して)最小ノルム解を与えるという特徴づけがあって、両者は、次元が無限か有限かということを除けば、本質的に同じものであるように見えた。のだけど、両者は、そこらに落ちてる説明を読んでも、見掛け上同じ計算をしているようには見えないので、ちゃんと同じに見えるようにしようと思ったのが、こういうことを考えた動機。なので、別に何か応用があるというわけではない。
まぁ、分かってみると当たり前のことのようではあるけど、要するに対称性がある特殊な"基底"(通常一次独立でないので、生成系とか呼ぶ方が正しい気がするけど)を使う場合には、擬似逆行列(最小ノルム解)を求める時に計算量が大きく減らせる(上に、多分計算精度も高い)というのがGabor変換etc.の大きな意義だということが分かる。一方で、工学的な問題で、そんな対称性の高い生成系を用いるのは、本当に正しいことだろうかという気もする
エルマンネットワーク
というのを実装して見た。リカレントニューラルネットワーク(RNN)の一形態？
単純再帰型ネットワーク
http://www.cis.twcu.ac.jp/~asakawa/waseda2002/elman.pdf
import numpy as np
def sigmoid(x):
return np.tanh(x)
def dsigmoid(x):
return 1.0-x*x
class ElmanNet:
def __init__(self , numIn , numHidden , numOut):
self.input = np.ones( numIn )
self.output = np.ones( numOut )
self.hidden = np.ones( numHidden )
self.context = np.zeros( numHidden )
self.weight_in_hid = np.random.uniform(low=-0.5,high=0.5,size=(numIn , numHidden) )
self.weight_hid_out = np.random.uniform(low=-0.5,high=0.5,size=(numHidden,numOut) )
self.weight_cont_hid = np.random.uniform(low=-0.5,high=0.5,size=(numHidden,numHidden) )
def propagate_forward(self , data):
self.input = np.copy(data)
self.context = np.copy(self.hidden)
self.hidden = np.vectorize(sigmoid)(np.dot(self.input , self.weight_in_hid) + np.dot(self.context , self.weight_cont_hid))
self.output = np.vectorize(sigmoid)(np.dot(self.hidden , self.weight_hid_out))
return self.output
def propagate_backward(self , target , learnRate=0.1):
error_out = (target - self.output) * np.vectorize(dsigmoid)(self.output)
error_hidden = np.dot(error_out , self.weight_hid_out.T) * np.vectorize(dsigmoid)(self.hidden)
self.weight_hid_out += np.array( learnRate * np.matrix(self.hidden).T * np.matrix(error_out) )
self.weight_in_hid += np.array( learnRate * np.matrix(self.input).T * np.matrix(error_hidden) )
self.weight_cont_hid += np.array( learnRate * np.matrix(self.context).T * np.matrix(error_hidden) )
def train(self , samples , r=1500):
L = len(samples)
for _ in xrange(r):
for (inp,outp) in samples:
self.propagate_forward(inp)
self.propagate_backward(outp)
def predict(self , data):
return self.propagate_forward(data)
if __name__=="__main__":
def tobit(n,nbits=6):
return [int(n==i) for i in xrange(nbits)]
net = ElmanNet(6,15,6)
samples = [(tobit(0),tobit(1)),(tobit(1),tobit(2)),(tobit(2),tobit(3)),(tobit(3),tobit(4)),
(tobit(4),tobit(5)),(tobit(5),tobit(4)),(tobit(4),tobit(3)),(tobit(3),tobit(2)),
(tobit(2),tobit(1)),(tobit(1),tobit(0)),(tobit(0),tobit(1)),(tobit(1),tobit(2)),
(tobit(2),tobit(3)),(tobit(3),tobit(4)),(tobit(4),tobit(5)),(tobit(5),tobit(4))]
net.train(samples)
out = net.predict(tobit(4))
print(out)
out = net.predict(out)
print(out)
print net.predict(tobit(2))
例として、0,1,2,3,4,5,4,3,2,1,0,1,2,3,4,5,...という時系列を学習させている。一応、値が上昇中か下降中か覚えておく必要があるけど、勝手に覚えておいてくれるという寸法。0,1,2,3,4,5は、[1,0,0,0,0,0],[0,1,0,0,0,0],[0,0,1,0,0,0],...にエンコードしてる。3bitで[0,0,0],[0,0,1],[0,1,],[0,1,1],...とかやった方がメモリには優しいけど、予測精度がでない
一応、上の計算結果は(乱数で初期化している部分があるので、結果は毎回異なるけど)例えば
[ 0.05445429 -0.0483922 -0.14895301 0.98319809 0.09161327 -0.09690022] [ 0.05653943 0.13345083 0.98843438 -0.28019387 -0.08493135 0.05609208] [-0.09952735 0.98482479 0.22479163 -0.13254797 -0.08240431 0.20785037]
とかなって、[0,0,0,1,0,0],[0,0,1,0,0,0],[0,1,0,0,0,0]に近い値が返ってくる
有限オートマトンとか、ある種の文法構造を学習することができるっぽいけど、他の文法推論アルゴリズムと比べて、大規模な文法の推論には向かないらしい
アルゴンの分子動力学計算
やったことなかったので、アルゴンの分子動力学計算をやってみた。まぁ、ただのGromacsの使い方メモのような。使用したGromacsのバージョンは、4.6.1
原理的なことに関しては、
分子動力学シミュレーションにおける温度・圧力制御第1回=能勢の熱浴と能勢・フーバー熱浴
https://www.jstage.jst.go.jp/article/mssj1998/10/44/10_44_29/_article/-char/ja/
分子動力学シミュレーションにおける温度・圧力制御 第3回：速度スケーリング法，ガウス束縛法，ベレンゼン熱浴
https://www.jstage.jst.go.jp/article/mssj/11/2/11_2_2_43/_article/-char/ja/
あたりを見ると、大体書いてある。まぁ、MDは原理自体は難しいことは、あまりないと思う。
ポテンシャルは、Lennard-Jonesポテンシャルで、あまり関係ないけど
量子力学はLennard-Jonesポテンシャルを再現できるか
http://d.hatena.ne.jp/m-a-o/20121217#p2
も参照。きっと、2030年くらいになれば、怪しい近似ポテンシャルとか力場とか使わなくても、量子力学的に原子核間に働くヘルマン・ファインマン力を計算して、原子核を古典論的に時間発展させるくらいのことは、そこらのPCで出来るようになるかもしれないけど、現在のところ厳しそう
(1)トポロジーファイルというものを作る。今の場合、実質的には、Lennard-Jonesパラメータと原子数の指定くらいしかやってない
;; argon.top [ defaults ] ; nbfuncは非結合粒子間ポテンシャルで、LJポテンシャルかBuckinghamポテンシャルを指定する(LJなら、1) ; comb-ruleは、異原子間のLJパラメータとかの決め方(sigma/epsilonを指定する場合は、2か3)。今はどっちでも同じ ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 3 no 1.0 1.0 [ atomtypes ] ; name mass charge ptype sigma(nm) epsilon(kJ/mol) AR 39.948 0.0 A 0.34 0.99774 [ moleculetype ] ; name nrexcl Ar 1 [ atoms ] ; id type resnr residu atom cgnr charge 1 AR 1 Ar Ar 1 0 [ system ] Argon [ molecules ] AR 512
(2)MD計算の初期状態を指定する。
pdbファイルで指定できるけど、その場合は、初速度の指定はできない(後で、Gromacsが適当な速度を付加してくれる)
速度も指定したい場合は、.groファイルフォーマットというのを使えばいいっぽい
とりあえず、何も考えず、立方格子の頂点に原子を配置するpdbを生成するコードをpythonで作る
#-*- coding:utf-8 -*-
#CellSizeの単位は(Å)
def gen_cubic_pdb(CellSize , N):
L = CellSize
A = CellSize/N
print('HEADER argon')
print('CRYST1 %9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d' % (L,L,L,90.0,90.0,90.0,"P 1",1))
atom_format = "ATOM %(serial)5d %(name)-4s %(resName)-4s%(chainID)1s%(resSeq)4d%(iCode)1s %(x)8.3f%(y)8.3f%(z)8.3f%(occupancy)6.2f%(tempFactor)6.2f"
for nx in xrange(N):
for ny in xrange(N):
for nz in xrange(N):
vars = dict()
vars['serial'] = nx*N*N+ny*N+nz+1
vars['x'] = nx * A
vars['y'] = ny * A
vars['z'] = nz * A
vars['name'] = 'Ar'
vars['resName'] = 'Ar'
vars['resSeq'] = nx*N*N+ny*N+nz+1
vars['iCode'] = ''
vars['occupancy'] = 1.0
vars['chainID'] = 'A'
vars['tempFactor'] = 0.0
print(atom_format % vars)
if __name__=="__main__":
gen_cubic_pdb(250.0 , 8)python genpdb.py > init.pdb
とかで、初期状態を指定する。トポロジーファイルに原子数を指定してあるので、一致させる
(3)MDの計算条件の指定(タイムステップやステップ数や相互作用のカットオフ等)
指定できる項目は、凄く多い
http://www.gromacs.org/Documentation/File_Formats/.mdp_File
とりあえず、最小に近い記述は、以下の通り(温度を273Kに保つという以外は、技術的な詳細)
;; NVT.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.005 ;; timestep(ps) 5fs nsteps = 200000 ;; total time=1000ps=1ns ; Output frequency and precision for xtc file nstxtcout = 1000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 273.0 ; Temperature coupling: no(NVE) , nose-hoover , berendsen tcoupl = berendsen tc-grps = system tau-t = 0.2 ;; time constant(ps) ref-t = 273.0 ;;reference temperature(K)
(4)MD計算
以上の入力ファイルを、まとめて、gromppというのにかけると、単一の入力ファイルを作ってくれる
grompp -p argon.top -c init.pdb -f NVT.mdp -o md.tpr
md.tprができれば成功。で、あとはmdrun
mdrun -deffnm md
出力は
md.xtc : trajectory
md.edr : エネルギー
md.gro : 最終的な位置と速度(継続して、別の計算に投げる時とかに使う)
md.log : ログファイル
など
(5)計算結果の確認と解析
g_energy -f md.edr -o energy -xvg none
とかやって、
End your selection with an empty line or a zero. ------------------------------------------------------------------- 1 LJ-(SR) 2 Disper.-corr. 3 Coulomb-(SR) 4 Potential 5 Kinetic-En. 6 Total-Energy 7 Temperature 8 Pres.-DC 9 Pressure 10 Vir-XX 11 Vir-XY 12 Vir-XZ 13 Vir-YX 14 Vir-YY 15 Vir-YZ 16 Vir-ZX 17 Vir-ZY 18 Vir-ZZ 19 Pres-XX 20 Pres-XY 21 Pres-XZ 22 Pres-YX 23 Pres-YY 24 Pres-YZ 25 Pres-ZX 26 Pres-ZY 27 Pres-ZZ 28 #Surf*SurfTen 29 T-System 30 Lamb-System
などと出てきたら、適当に
4 5 6 7 9
と打つと、energy.xvgというファイルが作られて、一定ステップごとに、ポテンシャルエネルギー・運動エネルギー・全エネルギー・温度・圧力の時間変化が見れる。とりあえず、温度が大体273Kに保たれていることとかが見て取れる。あと
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Potential -4.0542 0.18 1.72151 -0.99296 (kJ/mol) Kinetic En. 1739.87 0.031 1.07103 -0.172205 (kJ/mol) Total Energy 1735.81 0.2 1.36278 -1.1652 (kJ/mol) Temperature 273.003 0.0049 0.168055 -0.027022 (K) Pressure 1.23123 0.00017 0.00798183 -0.000493924 (bar)
とかいう平均値が出る。NVTの値から、理想気体の状態方程式によって、圧力を見積もると
p = NkT/V = 512*(1.3806488*1e-23)*273/(250*1e-10)^3 = 123508.202...(Pa)=1.23508...(bar)
くらいなので、大体合ってそう。Err.Est.が凄く小さいのは、多分最初からほぼ平衡状態にあることを示していると思う
MD計算で分かる物理量は、拡散係数、粘性係数(ずり粘性率、体積粘性率)、熱伝導率、比熱、密度、沸点、融点とか色々あると思うけど、一番簡単なのは、自己拡散係数で、平均二乗変位から見積もる方法と、速度相関関数からGreen-Kubo公式で求める方法があって、前者は、Gromacsに付いてるg_msdで簡単に計算できて
g_msd -s md.tpr -f md.xtc -o msd -xvg none
とかやって、適当に"0"と打つと、D=1.086e+04 (+/- 2985) 1e-5 cm^2/sであるとか返ってくる。
アルゴンガスの自己拡散係数のデータは探しても、あまり出てこないけど
DIFFUSION COEFFICIENT@thermopedia
http://www.thermopedia.com/content/696/
によると、常温常圧で、0.157(cm^2/s)とある。今、MD計算の方は、圧力を常圧にしてない(NPTアンサンブルでの計算は(6)を参照)けど、オーダーは合ってそう
g_msd -s md.tpr -f md.xtc -o msd -xvg none -beginfit 100 -endfit 600
とかやると、100psから600psのデータを使って自己拡散係数を出してくれるッぽい(今の場合、デフォルトでは100~900ps)。計算値は
100~400ps : 8346.9463 (+/- 2933.8926) 1e-5 cm^2/s 100~600ps : 9714.9385 (+/- 3332.4248) 1e-5 cm^2/s 100~800ps : 1.057e+04 (+/- 3164) 1e-5 cm^2/s 100~900ps : 1.086e+04 (+/- 2985) 1e-5 cm^2/s
1ns程度だと、まだ若干収束してない感じがする。計算時間を2nsと4nsに増やすと、それぞれ、1.242e+04 (+/- 1964) 1e-5 cm^2/sと1.287e+04 (+/- 166.5) 1e-5 cm^2/sくらいになって、標準偏差も大分小さくなるので、0.13(cm^2/s)弱くらいで収束するっぽい
Green-Kubo公式による自己拡散係数の計算もGromacsに付いているg_velaccとg_analyzeというので、できるらしい。
Diffusion Constant
http://www.gromacs.org/Documentation/How-tos/Diffusion_Constant
g_velaccというので、自己速度相関関数を計算してくれるっぽい。けど、そのためには中間速度を出力させる必要があって、それには、NVT.mdpに次の3行を追加して計算する。こうすると、md.trrとかいうファイルができるはず
nstxout = 100 nstvout = 100 nstfout = 100
計算できたら、
g_velacc -nonormalize -mol -f md.trr -s md.tpr g_analyze -f vac.xvg -integrate
とかで
Calculating the integral using the trapezium rule Integral 1 34.11745 +/- 0.00000
というのが出てくる。単位は、(nm^2/ps)らしいので、(cm^2/s)に直すためには、10^(-2)をかけて、これは速度自己相関関数の積分値なので、拡散係数にするために、3で割ると、0.1137(cm^2/s)という値になる。
(6)沸点・融点の計算
今度は、NPT一定の条件で、ちょっとずつ温度を下げていって、気液相転移・固液相転移を見る。歴史的には、Alderが、2次元系で剛体球ポテンシャルで、固液相転移を確認したのが、分子動力学の始まりらしい。
topファイルとpdbファイルは、さっきと同じで、mdpを以下のように変える
;; anneal.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.001 ;; timestep(ps) 1fs nsteps = 40000000 ;; total time=40ns ; Output frequency and precision for xtc file nstenergy = 10000 nstlog = 10000 nstxtcout = 10000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 150.0 ;; temperature and pressure coupling tcoupl = berendsen tc-grps = system tau_t = 0.2 ;; time constant(ps) ref_t = 150.0 ;;reference temperature(K) pcoupl = berendsen pcoupltype = isotropic nstpcouple = 1 tau_p = 0.015 ;; time constant(ps) compressibility = 4.5e-5 ;; (1/bar) ref_p = 1.0 ;;reference pressure(bar) ; SIMULATED ANNEALING annealing = single ; Type of annealing for each temperature group (no/single/periodic) annealing_npoints = 11 ; Number of time points annealing_time = 0 2800 3000 4800 5000 14800 15000 24800 25000 34800 35000 ; List of times(ps) at the annealing points annealing_temp = 150 150 110 110 90 90 70 70 50 50 30 ; Temp. at each annealing point
このように書くと、0~2800psでは、熱浴の温度を150Kに保ち、2800~3000psの間に、150Kから110Kへ冷却し(熱浴の温度をリニアに下げていく)、3000~4800psでは、温度を110Kに保ち...という感じの計算をやってくれる。圧力も1barに保つように設定されている。で、以下のコマンドで計算
grompp -p argon.top -c init.pdb -f anneal.mdp -o md.tpr mdrun -deffnm md
計算が終わると、今度はg_energyコマンドで、密度や体積も出力できるようになっている。各温度で平衡状態に於ける密度を見ると
150K:3.1(kg/m^3) 70K:7.4(kg/m^3) 50K:1690(kg/m^3) 30K:1742(kg/m^3)
くらいになっていて、70~50Kの間で、密度の急激な増大が起きていることが確認できる。実験値は
気体アルゴン密度(0℃、1気圧):1.784(kg/m^3) 液体アルゴン密度(87K、液体):1380(kg/m^3) 固体アルゴン密度(40K、固体):1656(kg/m^3) 融点:84K 沸点:87K
なので、融点・沸点は全然ずれている。密度は、大体合っているっぽい。一応、最終配置から、30Kに保って、自己拡散係数を出してみると、液体でなく、固体になっているっぽい
逆に、固体から始めて温度をあげていく。以下のようなpythonコードで、fcc結晶のpdbを作成する。今度は、原子数は3x3x3x4=108個
#!/usr/bin/env python
def gen_fcc(a , Nx , Ny , Nz):
ret = set([])
for n1 in xrange(2*Nx):
for n2 in xrange(2*Ny):
for n3 in xrange(2*Nz):
if (n1+n2+n3)%2==1:continue
x = n1*0.5*a
y = n2*0.5*a
z = n3*0.5*a
yield (x,y,z)
if __name__=="__main__":
a = 5.26 #-- lattice constant
Nx,Ny,Nz = 3,3,3
print('HEADER argon')
print('CRYST1 %9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d' % (Nx*a,Ny*a,Nz*a,90.0,90.0,90.0,"P 1",1))
atom_format = "ATOM %(serial)5d %(name)-4s %(resName)-4s%(chainID)1s%(resSeq)4d%(iCode)1s %(x)8.3f%(y)8.3f%(z)8.3f%(occupancy)6.2f%(tempFactor)6.2f"
for n,(x,y,z) in enumerate(gen_fcc(a,Nx,Ny,Nz)):
vars = dict()
vars['serial'] = n+1
vars['x'] = x
vars['y'] = y
vars['z'] = z
vars['name'] = 'Ar'
vars['resName'] = 'Ar'
vars['resSeq'] = n+1
vars['iCode'] = ''
vars['occupancy'] = 1.0
vars['chainID'] = 'A'
vars['tempFactor'] = 0.0
print(atom_format % vars).mdpは以下のように適当に作った(.topファイルの原子数も108に変更しておくこと)
;; heatup.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.001 ;; timestep(ps) 1fs nsteps = 40000000 ;; total time=40ns ; Output frequency and precision for xtc file nstenergy = 100000 nstlog = 100000 nstxtcout = 100000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency rlist = 0.5 ;;(nm) rcoulomb = 0.5 rvdw = 0.5 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 50.0 ;; temperature and pressure coupling tcoupl = berendsen tc-grps = system tau_t = 0.2 ;; time constant(ps) ref_t = 50.0 ;;reference temperature(K) pcoupl = berendsen pcoupltype = isotropic nstpcouple = 1 tau_p = 0.015 ;; time constant(ps) compressibility = 4.5e-5 ;; (1/bar) ref_p = 1.0 ;;reference pressure(bar) ; SIMULATED ANNEALING annealing = single ; Type of annealing for each temperature group (no/single/periodic) annealing_npoints = 11 ; Number of time points annealing_time = 0 2800 3000 4800 5000 14800 15000 24800 25000 34800 35000 ; List of times(ps) at the annealing points annealing_temp = 50 50 60 60 80 80 100 100 120 120 160 ; Temp. at each annealing point
今度は、平衡状態に達した後で、各温度での密度が
50K:1750(kg/m^3) 80K:1680(kg/m^3) 100K:1620(kg/m^3) 120K:1290(kg/m^3) 160K:3.01(kg/m^3)
くらいになっているので、100~120Kで液体化して、120~160Kの間で気体化しているっぽい。厳密には、もっと細かく点を取るべきだと思うけど、冷却していった時と、加熱していった時で、沸点と融点が全然合ってそうにないのは、何故だろう？一応、これらの値は、原子数を500に増やしても、ほぼ同一だった
(7)疑問
心当たりのある方は、教えてください
・沸点や融点が実験値とかなりずれる理由は何？(単純に、ポテンシャルのfittingがよくない?)
・固体から加熱していった時と、気体から冷却していった時で、融点・沸点が大きくずれる理由は何？経験的に知られているように、固体の液化温度と液体の固化温度が概ね一致するための条件は？
