AbInitioMC
平行移動(shift)と回転(rotation)を抑える
モンテカルロ試行を実行している最中に系全体の平行移動(shift)と回転(rotation)が起こる． これを抑える必要がある．これによって止められる自由度は系全体の自由度に対してたかがしれているので，十分に自由エネルギーを計算できると考えられる．
回転は各座標でのモーメントによって計算できる．
当初は系全体のshiftとrotationの自由度6を抑えるために，二つの原子を選んで，解析的に解こうと考えた．しかし，この実装においていくつかの問題が出てきて，その過程で方針を変更した．
- moment計算の中心をどこに取るか？
- momentとshiftを解析的に解いた場合に，連立方程式の係数行列が反対称となり，解がうまく求まらない．
- そもそもshiftしてしまうとmomentが変わる．
- 従って，解は2回に分けて求める必要がある．
- 最初に，中心を戻す
- momentを消す．
- ならば，系全体を回した方が楽かも．
という事で，系全体を戻すこととした．なお，2原子による制約の試行錯誤の様子は，TwoAtomsSuppressingにまとめている．
rotationの特徴
ベクトルのrotationは反対称行列を生成する．これは，正則ではなく逆行列を持たない．したがって，3つのベクトルからもとまるキャンセル操作を行う必要がある．現在のところ，
initial.rbによって移動させてみる
[bob-no-MacBook-Pro:~/Desktop/FP-MD&MC/nvt_mc_sample_2] bob% ruby initial.rb 移動距離の上限を，原子間距離に対する比で入力してください． 0.02 "POSCAR" Center position of the system: 0.5004517388 0.5004170598 0.5000285272 "Al_POSCAR" Total moments of the system: -0.0024100926 0.0005558295 0.0080836649 7 18 15 [0.375, 0.375, 0.125] [-0.000577614999999976, 0.003299263, 0.001716415] [0.375, 0.625, 0.875] [0.000939715000000008, -0.00195212099999997, -0.00241799099999995] [0.375, 0.125, 0.375] [0.00216666799999998, 0.001630883, 0.00029815900000002] rxf:=Vector([ -0.0024100926, 0.0005558295, 0.0080836649]): p1:=Vector([ 0.3750000000, 0.3750000000, 0.1250000000]): f10:=Vector([ -0.0005776150, 0.0032992630, 0.0017164150]): f1:=Vector([dx, 0.0032992630, 0.0017164150]): p2:=Vector([ 0.3750000000, 0.6250000000, 0.8750000000]): f20:=Vector([ 0.0009397150, -0.0019521210, -0.0024179910]): f2:=Vector([ 0.0009397150,dy, -0.0024179910]): p3:=Vector([ 0.3750000000, 0.1250000000, 0.3750000000]): f30:=Vector([ 0.0021666680, 0.0016308830, 0.0002981590]): f3:=Vector([ 0.0021666680, 0.0016308830,dz]):
Moment_cancel4.mwで計算
restart;
cc:=Vector([0.5,0.5,0.5]):
rxf:=Vector([ 0.0004089012, 0.0000018379, -0.0001244011]):
p1:=Vector([ 0.8750000000, 0.1250000000, 0.3750000000]):
f10:=Vector([ -0.0000053570, -0.0001813910, -0.0000237270]):
f1:=Vector([dx, -0.0001813910, -0.0000237270]):
p2:=Vector([ 0.1250000000, 0.8750000000, 0.8750000000]):
f20:=Vector([ 0.0001444470, -0.0000486220, -0.0000577100]):
f2:=Vector([ 0.0001444470,dy, -0.0000577100]):
p3:=Vector([ 0.1250000000, 0.1250000000, 0.1250000000]):
f30:=Vector([ 0.0000235240, 0.0001876580, 0.0000114410]):
f3:=Vector([ 0.0000235240, 0.0001876580,dz]):
with(LinearAlgebra):
dd:=Vector([dx,dy,dz]);
Vector[column](%id = 18446744078172353950)
pf1:=CrossProduct(p1-cc,f1)-CrossProduct(p1-cc,f10);
pf2:=CrossProduct(p2-cc,f2)-CrossProduct(p2-cc,f20);
pf3:=CrossProduct(p3-cc,f3)-CrossProduct(p3-cc,f30);
Vector[column](%id = 18446744078172355870)
Vector[column](%id = 18446744078184022974)
Vector[column](%id = 18446744078184024654)
eq1:={seq(pf1[i]+pf2[i]+pf3[i]+rxf[i]=0,i=1..3)};
{-HFloat(1.4062547500000002e-4) + HFloat(0.375) dx
- HFloat(0.375) dy = 0,
HFloat(3.94958325e-4) - HFloat(0.375) dy - HFloat(0.375) dz = 0,
-HFloat(0.125) dx - HFloat(3.1220999999999993e-6)
+ HFloat(0.375) dz = 0}
s1:=solve(eq1,{dx,dy,dz});
{dx = 0.001064923400, dy = 0.0006899221333, dz = 0.0003633000667}
sx:=dx=subs(s1,dx);
tmp:=p1+subs(sx,f1);
printf("%15.10f %15.10f %15.10f\n",tmp[1],tmp[2],tmp[3]);
dx = 0.001064923400
Vector[column](%id = 18446744078184025854)
0.8760649234 0.1248186090 0.3749762730
sy:=dy=subs(s1,dy);
tmp:=p2+subs(sy,f2);
printf("%15.10f %15.10f %15.10f\n",tmp[1],tmp[2],tmp[3]);
dy = 0.0006899221333
Vector[column](%id = 18446744078184018878)
0.1251444470 0.8756899221 0.8749422900
sz:=dz=subs(s1,dz);
tmp:=p3+subs(sz,f3);
printf("%15.10f %15.10f %15.10f\n",tmp[1],tmp[2],tmp[3]);
dz = 0.0003633000667
Vector[column](%id = 18446744078184019478)
0.1250235240 0.1251876580 0.1253633001
POSCARへいれて修正
analysisの結果は
[bob-no-MacBook-Pro:~/Desktop/FP-MD&MC/nvt_mc_sample_2] bob% ruby analysis.rb POSCAR Al_POSCAR "POSCAR" Center position of the system: 0.4999310388 0.5019172761 0.4983274699 "Al_POSCAR" Total moments of the system: -0.0000000000 -0.0000000000 -0.0000000000
このときPOSCARはこんなかんじ．
New structure
1.0
8.0827999115 0.0000000000 0.0000000000
0.0000000000 8.0827999115 0.0000000000
0.0000000000 0.0000000000 8.0827999115
32
Direct
0.121522145 0.123875609 0.126380623
0.128159579 0.126352422 0.626918472
0.127259946 0.624477652 0.127559410
0.124400544 0.626117690 0.628206277
0.628719727 0.124872134 0.125721671
0.627014433 0.127880424 0.627861575
0.628328788 0.625543923 0.122662104
0.3577599858 0.3782992630 0.1267164150
0.372341665 0.371619974 0.626845366
0.378018697 0.878015699 0.125673938
0.374832130 0.877344816 0.624057103
0.874755976 0.377612707 0.122731063
0.876312461 0.371798808 0.626256769
0.874549299 0.873497035 0.121960413
0.871449574 0.878544546 0.624740749
0.3771666680 0.1266308830 0.3208643254
0.378196507 0.123765858 0.877533227
0.371450778 0.628124147 0.373329702
0.3759397150 0.6710547990 0.8725820090
0.873398263 0.127524271 0.377487908
0.877588955 0.122796382 0.871737531
0.877304117 0.626738136 0.377999761
0.872379884 0.627683811 0.871675476
0.128488889 0.371784581 0.372956765
0.124846778 0.375222287 0.874227586
0.127467735 0.871624886 0.376825300
0.122699096 0.877560089 0.873862379
0.628775212 0.373817284 0.372519349
0.622768518 0.372981206 0.875940198
0.625368467 0.878075847 0.371577766
0.624286663 0.872248360 0.872548139
0.624242046 0.627867306 0.628519669
Keyword(s):
References: