2012-02-23

剛体球 MD

剛体球で分子動力学シュミレーションやった。
衝突はめり込みで判定して、セルで切って次に衝突する粒子を隣接したセルからのみ探すようにしてる。

２次元10000粒子周期境界で計算。
小さくて見え辛いけど良く見ると結晶化してるのが見える。

こっちは完全に遊びの計算。たぶん8000粒子位。

2011-12-21

レナード・ジョーンズ

レナード・ジョーンズに従う粒子のシミュレーション。重力多体系のコードを流用して作った。
馬鹿正直に全ての粒子の相互作用を計算してるから無茶苦茶遅い。

[12/22追記]運動量の大きさで色分け。

質量の大きさを変化させることができるようにした。

2011-11-26

emacsでpdfファイルを閲覧する

emacs

latexでのpdfの作成〜閲覧までをemacsの中で行いたかったので調べてみました。

doc-viewを使えばcocoa emacsでpdfを表示させることが出来るようです。
Emacs 23でPDFを表示させる--新機能「DocView」を試す - builder by ZDNet Japan

少し設定が必要見たいですが、僕の環境ではとくに何もしないでもemacsでpdfファイルを開くとdoc-viewが起動しました(!?)

latexで編集しながらpdfファイルを閲覧する際には ~/.emacsに

(add-hook 'doc-view-mode-hook 'auto-revert-mode)

と書いて自動更新をonにしておくと便利です。

編集の様子はこんな感じです。pdfの更新とファイルのリロードに数秒のタイムラグがあるのが若干気になります、、、

2011-11-06

yatex-modeになったときにマイナーモードcdlatexを起動する

emacs latex

yatexそのものにも補完機能はあるのですが普段org-mode+cdlatexを使っているでそちらの方に手が慣れいます。
そういうわけで、yatex起動時にcdlatexも起動する設定。

(require 'cdlatex)
(add-hook 'yatex-mode-hook 'turn-on-cdlatex)

そういえば、cdlatexってもともと入ってんだっけ?

2011-10-30

Org-modeでBibtexを使う

latex org-mode emacs

Org-modeからbibtexを使う手順の覚え書き。
BibTeX と Org Export | Amrtaを参考にしました。

bibtexファイルをカレンドディレクトリに放り込んで
orgファイルの先頭に

#+BIBLIOGRAPHY bibtex file name

と書き、文献を表示したい場所(ファイル末とか)に

\bibliographystyle{plain}
\bibliography{bib tex file name}

と記入すれば良いようです。

後は引用したい場所で、C-cC-xからキーワードを入力すれば文献と対応した番号が表示されます。

biographyを二回書かないといけないのが面倒で多分ちゃんとしたやり方があるんだろうけど、とりあえず動いたので
いいや。

2011-10-19

負のエントロピー?

突然ですが理想気体のエントロピーを求めてみましょう。

理想気体分配関数は、
$Z=\frac{1}{h^{3N}}d^{3N}pd^{3N}p\exp(-\beta \frac{p^2}{2m})$
$=V^N(\frac{2\pi m}{h^2\beta})^{3N/2}$

自由エネルギーは、 $F=-\log(Z)/\beta$ これからエントロピーは、
$S=\frac{\partial F}{\partial T} = k_B(\log(V(\frac{2\pi m k_B T}{h^2})^{3/2})-1)$

この表式を見ると温度を下げていくとエントロピーが負になってしまうように見えます。
また、T=0では $S\to \infty$ となってしまいます。

一見不合理なこの結果は、不確定性関係を忘れてしまったことによります。
logの中身を見てみると、分子に $V m k_B T \approx \Delta x^3\Delta p^3$ 分母に hが出て来ています。
不確定性関係 $\Delta p\Delta x > h$ の制限を考慮すれば、単純に $T\to 0$ としたからといってlogの中身を0にすることは出来ないわけです。

これは分配関数を計算する際に出てきた係数 $h^{3N}$ が位相空間での状態1つに対応する体積だったことからも理解できます。

エントロピーが負になるという現象は確率変数が連続的な値を取り得る場合に見られて例えばDifferential entropy - Wikipediaに例が載っています。

2011-10-05

揺動散逸定理についてつらつらと

物理統計力学 matplotlib python

少し息抜きをかねて揺動散逸定理について感じたことをつらつらと。
揺動散逸定理というのは、アインシュタインが発見した定理で、
平衡状態から離れるよう働く揺動と平衡状態へと向かわせる散逸の間の関係を表したものです。

式で表せば、
$\frac{dv}{dt} = -\frac{\partial H}{\partial x}-\mu v(t) + \delta F(t)$
$<\delta F(t) \delta F(t')> = 2D\delta(t-t')$
この、　散逸の大きさ $\mu$ と、　揺動力の大きさ $D$ の間に、 $D=\mu k_B T$ の関係が成立するというものです。

導出は幾通りかあって、
ひとつは、系が最初から平衡状態にあるとして速度分散( $k_BT/m$ )が時間に依らないためには上の揺動散逸定理が成り立たなければならないというもの。

もうひとつは、最初は系が平衡状態にないとするものです。
十分時間が経てば速度分散は揺動と散逸の釣り合いによって初期条件に依らない一定値に落ち着くだろうと予想されます。この時の速度分散がある温度での平衡状態での速度分散と等しいと置くことによって揺動散逸定理を導出します。

二つ導出方法は、時間の原点を注目する粒子が平衡状態に達する前に置くか後に置くかという違いを無視すれば実質的に同じものです。

ここでは、二つ目の導出について簡単なシミュレーションを行ってみようと思います。

調和振動子の場合を考えます。　揺動と散逸を固定して適当な初期条件を与えて位相空間上で時間発展させたものをプロットしてみます。

from pylab import *
time = range(10000)
for i in range(10):
    del(x,v)
    x = [2.0*(random()-0.5)]
    v = [2.0*(random()-0.5)]
    dt = 0.01

    for t in time:
        x += [x[-1] + v[-1] * dt]
        v += [v[-1] - x[-2] * dt
              - v[-1] * dt
              + 5*(random()-0.5) * dt]
    plot(x, v)

１０本の軌跡が図の中心に引き寄せられて最終的には初期条件に依らない一定軌道を描きます。
この軌道上での速度分散が、揺動と散逸がない場合のハミルトニアンのある温度での速度分散と等しいというのが揺動散逸定理の主張するところです。

揺動と散逸が熱浴の役割を果たしているということもできるでしょう。

大分簡単化されたモデルではありますが、初めて非平衡状態から平衡状態への緩和の過程を記述したという点でこの定理のもつ意義は非常に大きかったようです。