|
|
> No.4247[元記事へ]
SECONDさんへのお返事です。
（５）NVEの場合
> それよりもっと、大事な事がありました。
>
> 初期の温度から、外れたところで安定する原因について。
>
> 開始の particle 配置による、ポテンシャル・エネルギー総和 と、
> 運動中の particle 位置による、ポテンシャル・エネルギー総和の平均 との差分が
> 運動エネルギーとして、particle の初期温度に、加算されるようです。
> 開始の particle 配置による、ポテンシャル・エネルギー総和 と、
> 運動中の particle 位置による、ポテンシャル・エネルギー総和の平均 との差分が
> 運動エネルギーとして、particle の初期温度に、加算されるようです。
これは、本プログラムのような断熱系(NVE:粒子数一定、体積一定、全エネルギ一定）のシステムでは
正常な振る舞いです。
分子間のポテンシャルはsigmaより少し離れるとV(r) ～ r^(-6)で浅くなりますので、
初期状態を等間隔にするとどうしても平均のポテンシャルエネルギーは大きくなり、
粒子が動いて近くに来ることが多くなると、運動エネルギーに変わって、温度が上がります。
簡易的に温度を一定にする方法として、時々ajustVelocity(temp)を呼ぶという方法もあります。
SUB ajustVelocity(temp)
LET w=SQR( temp / sysTemp)
MAT vxy= w*vxy
MAT vz= w*vz
END SUB
（６）温度を返す関数として
FUNCTION sysTemp
LET v2 = 0
FOR i=1 TO Nmt
LET v2 = v2 + ABS(vxy(i))^2-vz(i)^2 ! vx^2 + vy^2 + vz^2
!LET v2 = v2 + ABS(ABS(vxy(i))+vz(i))^2 ! vx^2 + vy^2 + vz^2
NEXT i
LET sysTemp = 0.5*mass*v2/Nmt/k_ !(totalEnergy/Nmt)/k
END FUNCTION
を定義されていますが、
統計力学によると、粒子1個あたり、1自由度あたり (1/2)*k*TTempのエネルギーを分配されますから
2D : Ek = (2/2)k*Temp*Nmt ,
3D : Ek = (3/2)k*Temp*Nmt
となりますので、3Dの場合は Temp = (2/3)*Ek/(Nmt*k) となります。
|
|