我們分別用實驗和模擬的方法觀察非平衡穩定態系統，用高速攝影機追蹤在 振動台提供垂直振動下被氣體顆粒環繞的顆粒布朗物體，觀察其行為。發現其速 度機率分布為高斯分布，並且遵守能量均分定律。接下來我們希望在模擬中找出 能使顆粒布朗物體滿足能量均分的情況，模擬時我們觀察二維系統，並且用溫度 調整的方式代替振動台的作用。在實驗中我們可以量測垂直碰撞係數，和摩擦係 數，將這些參數代入分子動力學模擬中，但是我們發現無論是顆粒氣體或是顆粒 布朗物體都不遵守能量均分定律，之後我們改變摩擦係數和垂直碰撞係數，發現 除非摩擦係數 大於 50，否則粒子不會遵守能量均分定律，它每自由度的平均轉 動動能都低於平均移動動能。
我們的模擬跟前人作的模擬的差別是我們的顆粒布朗物體，它的質量、半徑 都大於顆粒氣體，並且它沒有獲得任何的外加能量。相對而言，顆粒氣體有獲得 外加的能量。這樣考慮的理由是因為我們假設在實驗中顆粒布朗物體的水平運動 幾乎不受到振動台的影響，它的水平速度的產生主要是經由跟顆粒氣體碰撞，所 以在模擬中對於整個顆粒體系統而言，顆粒布朗物體是一個特例。在這種情況下，
我們希望能得到跟前人作的模擬不相同的結果，也就是我們實驗中得到的顆粒布
朗物體滿足能量均分定律之結果。不過最後的結果似乎跟大部分文獻相同，能量
是不均分的。而且模擬中在摩擦係數小的情況下，顆粒布朗物體的轉動動能是很
33
低的，也就是顆粒布朗物體幾乎不會轉動。但是明顯違反我們的實驗數據。事實 上在實驗過程中，光用肉眼就可以看到顆粒布朗物體的轉動是很劇烈的，因此我 們的模擬明顯沒有辦法完全代表實驗的情況。
而這篇論文中實驗跟模擬最大的不同點是實驗的情況是三維，但我們模擬只
模擬二維，三維時顆粒布朗物體被碰撞機率將大幅提升。另外在實驗中顆粒布朗
物體下方的顆粒氣體，在我們的模擬之中是沒有考慮的。雖然顆粒布朗物體並未
直接受到振動台的作用，但或許間接的受到其下方的顆粒氣體撞擊而產生水平運
動和轉動，如果我們能確認使顆粒布朗物體產生的水平速度和轉動的主要機制為
何，就可以解釋我們實驗結果和模擬結果的不同點。
34
圖 2-1
圖 2-1 實驗儀器圖:
圖(a)本實驗所用的振動台為Labworks 公司的型號為(LW 127 123-500)振動台機組，
最大振動頻率可達 4500Hz。圖(b)為用來攝影的高速攝影機的型號 Prosilica GE680 CCD camera 和 IDT MotionPro X3。
(a)
(b) 33 公分
35
圖 2-2 顆粒布朗物體結構示意圖:
(a)顆粒布朗物體的製成示意圖，將壓克力材料挖空後後，將磁鐵裝入，最後將有黑 點的白紙貼上，最後在側邊貼上黑色膠帶，完成圖如(b)。
圖 2-2
有黑點的白紙
壓克力材料
黑色膠帶 (a)
(b)
強力磁鐵
5 公分 (b)
36
高速攝影機
容器
顆粒布朗物體
顆粒氣體
圖 2-3 實驗儀器架設圖
(a)、(c)是實驗儀器照片，將容器(紅框部分)固定在振動台上，並將高速攝影機架 設於上面。(b)、(d)為示意圖，容器中裝入顆粒布朗物體(有兩黑點的白色圓形部分) 和顆粒氣體(黃色圓形部分)、並在容器周圍加了一圈方形強力磁鐵，北極朝上。
(c) 圖 2-3
(a) (b)
(d)
高速攝影機
振動台
50cm 5cm
強力磁鐵 強力磁鐵 容器
37
圖 2-4 轉動慣量測量圖:
以圓形為例，圖中灰色圓盤為顆粒布朗物體，懸吊線(紅線)的長度為 L ，假設顆粒布
朗物體的質量為m，物體的外接圓半徑為r，我們使顆粒布朗物體做小角度扭動，並
記錄其轉動週期
，就可以根據右上公式計算顆粒布朗物體的轉動慣量 I ，右方的表格為對於測量的方法和積分的方法的結果比較，積分的時候假設顆粒布朗物體是均勻 的，不考慮中間挖洞裝入強力磁鐵對於轉動慣量所造成的誤差。
圖 2-4
形狀 測量的轉動慣量
(克‧平方公分)
積分的轉動慣量 (克‧平方公分)
三角形 30.84.4 36.7
圓形 155.79.3 155.6
正方形 122.28.3 121.8
正六邊形 97.96.6 109.8
𝐼 = 𝑚𝑔𝜏 2 𝑟 2 4𝜋 2 𝐿
顆粒布朗物體
2𝑟 𝐿
38
圖 2-5 碰撞係數、摩擦係數量測示意圖:
(a)為碰撞係數量測方法，黃色的圓形為顆粒氣體，讓顆粒氣體從高度
h
0落下後，並落到藍色的基板上並反彈，假設虛線為顆粒氣體落下之軌跡(只考慮垂直方向)，第
i
次反彈後高度
h
i，可用圖中公式求得顆粒氣體與基板的垂直碰撞係數。(b)摩擦力量測示意圖，增加傾斜角度，讓兩板剛好開始產生相對運動，記錄此時的角度，顆
粒氣體和顆粒布朗物體的摩擦係數的摩擦係數
1以及顆粒氣體彼此之間
2可以表示為
𝜇
1,2= 𝑡𝑎𝑛𝜃
。基板
𝜇
1,2= 𝑡𝑎𝑛𝜃
(a)
(b)
圖 2-5
顆粒氣體
底板
𝛼
1,2= √ ℎ
iℎ
𝑖−139
圖 2-6 顆粒布朗物體影像分析圖:
圖(a)為顆粒布朗物體的照片，顆粒布朗物體在半徑為 25 公分的圓型容器中運動，如 圖 2-3(c)。在運動過程中我們追蹤它的位置及角度。圖(b)為影像分析示意圖，顆粒 布朗物體上的白紙畫的兩個黑點 A、 B ，經影像處理後，顆粒布朗物體的白色部分中
心即顆粒布朗物體的x、 y 座標， A 、 B 兩黑點可以決定顆粒布朗物體的方向(或角
度)。
圖 2-6
( , )
B
θ
A
顆粒布朗物體
(a) (b)
40
圖 2-7 計算顆粒布朗物體速度示意圖:
用高速攝影機以 100 張/秒的速度連續拍攝 10 張影像，觀察顆粒布朗物體在 0.1 秒
內的(a)x位置、(b) y 位置、(c)角度隨時間變化，時間間隔為 0.01 秒，此時間間隔
中顆粒布朗物體的速度在彈道特徵，我們求出擬合直線，其斜率即為顆粒布朗物體的 速度。
圖 2-7
時間(秒)
時間(秒)
時間(秒) 位
置
(
公 分)
位 置
(
公 分)
角 度
(
弧 度)
(a)
(b)
(c)
41
圖 2-8
圖 2-8 顆粒氣體位置、速度分析圖:
用高速攝影機以 500 張/秒的速度連續拍攝 5 張在 1.42x1.42 平方公分區間中的顆粒氣體的 照片如(a)，追蹤其位置中心((a)圖紅色十字)，並留下良好的數據。得到 5 張影像中的顆粒 氣體位置軌跡如圖(b)，紅色的十字點為顆粒氣體中心位置，把(b)圖中籃色虛線圈起的顆粒
氣體的x、 y 位置分別對時間作圖，得到(c)、(d) 圖，算出(c)、(d)圖中擬合直線的斜率即
為顆粒氣體的x、 y 方向速度，分別為-7.91 公分/秒、-2.37 公分/秒。
(a) (b)
(c) (d)
位 置
(
公 分)
時間(秒) 時間(秒)
位 置
(
公 分)
位置
(
公 分)
x位置(公分)
42
圖 2-9 顆粒布朗物體位置分布圖:
我們將前 100 秒顆粒布朗物體的位置標註在圖中，圖中藍色圓圈為容器的邊界，可 以發現顆粒布朗物體因受磁力作用分布在一定範圍。
y 位 置
(
公 分)
圖 2-9
(b)
x位置(公分)
43
圖 2-10 顆粒布朗物體位置，角度隨時間變化圖:
(a)100 秒內每隔 0.5 秒顆粒布朗物體的x、 y 座標隨著時間變化圖，(b)100 秒內每
隔 0.5 秒顆粒布朗物體角度隨時間變化圖。
圖 2-10
(a)
角 (b) 度
(
弧 度)
時間(秒)
時間(秒)
𝒙
𝒚位 置
(
公 分)
44
2-11 顆粒布朗物體速度隨時間變化圖:
100 秒內(a)顆粒布朗物體每隔 0.5 秒x、 y 方向速度隨時間的變化圖。(b)顆粒布朗
物體每隔 0.5 秒角度隨時間的變化圖。
圖 2-11
(a)
(b) (b) 速
度
(
公 分/
秒)
時間(秒)
時間(秒) 角
速 度
(
弧 度/
秒)
𝑽
𝒙 𝑽𝒚45
圖 2-12 顆粒布朗物體的能量變化:
顆粒布朗物體的三個自由度的能量在 100 秒內每隔 0.5 秒隨時間變化圖，其中
E
x是x方向動能，
E
y是y
方向動能、E
是轉動動能。圖 2-12
動 能
(
微 焦 耳)
時間(秒) 𝑬𝒙
𝑬𝒚 𝑬𝝎
46
圖 2-13 顆粒布朗物體能量自相關圖:
6 秒內的能量自相關函數圖，橫軸為時間，黑色線為x方向動能自相關函數，紅色
線為 y 方向動能自相關函數，藍色線為轉動動能的自相關函數。
圖 2-13
能 量 自 相 關 函 數
時間(秒)
方向 方向 轉動
47
圖 2-14 顆粒布朗物體、顆粒氣體的速度機率分布圖:
(a)顆粒布朗物體速度的機率分布圖，Vx是x方向速度，Vy是
y
方向速度，單位是(公分/秒)，是角速度，單位是(弧度/秒)，(b)顆粒氣體移動方向的速度機率分布，Vx
是x方向速度，Vy是
y
方向速度，實線為高斯分布函數擬合。(a)
(b)
圖 2-14
速度(公分/秒)、角速度(弧度/秒)
機 率
速度(公分/秒)
𝑽𝒙 𝑽𝒚 𝑽𝒙 𝑽𝒚 𝝎
48
圖 2-15 顆粒布朗物體和顆粒氣體能量分布機率圖:
(a)顆粒布朗物體(b)顆粒氣體之不同自由度之能量機率分布，
E
x是x方向動能，E
y是
y
方向動能、E
是轉動動能。圖 2-15
機 率
動能(微焦耳) 動能(微焦耳) (a)
(b)
機 率
𝑬𝒙 𝑬𝒚 𝑬𝝎
𝑬𝒙 𝑬𝒚
49
圖 2-16 不同振動強度
的圓形顆粒布朗物體平均能量圖:質量為 59.7 克、半徑為 2.5 公分、轉動慣量為 155.62 克‧平方公分的圓形顆粒布朗 物體在不同的振動強度下的三個自由度之平均動能。Ex 是x方向平均動能，Ey 是
y
方向平均動能、E 是轉動平均動能。𝑬𝒙 𝑬𝒚 𝑬𝝎
圖 2-16
振動強度Г
動 能
(
微 焦 耳)
50
圖 2-17 不同形狀、質量顆粒布朗物體平均能量:
不同形狀、質量的顆粒布朗物體在相同的振動台振動強度
值下，三個自由度之每單位質量的平均動能。Ex 是x方向平均動能、Ey 是
y
方向平均動能、E 是轉 動平均動能。𝑬𝒙 𝑬𝒚 𝑬𝝎
圖 2-17
質量(克) 單
位 質 量 平 均 動 能
(
微 焦 耳/
克)
51
表 2-1 圓盤顆粒布朗物體參數和實驗結果:
圓形顆粒布朗體的參數，以及振動台的振動強度振動強度，和實驗完成後的平均動 能。
表 2-1
參數
物體質量 59.7(克)
物體底面半徑 2.51(公分)
物體高度 2.01(公分)
物體轉動慣量 155.62(克‧公分平方)
振動台振動強度 3.63
平均動能
E
x :x方向平均動能 67.3
2.0(微焦耳)
E
y :y方向平均動能 69.5
2.1(微焦耳)
E
:轉動動能 67.2
2.1(微焦耳)52
圖 3-1 分子動力學模擬能量自相關函數圖:
(a)能量自相關函數圖，橫軸為時間，黑色線為x方向移動動能，紅色線為 y 方向移
動動能，藍色線為轉動能量之自相關函數，(b)縱軸取對數後之圖形。
圖 3-1
時間(0.5 秒) 能
量 自 相 關 函 數
x方向 y方向
轉動
時間(0.5 秒) 能
量 自 相 關 函 數
(a)
(b) x方向
y方向 轉動
53
圖 3-2 顆粒體碰撞示意圖:
有旋轉的粗糙顆粒體 1,2 互相碰撞，黑點為球心，
𝑣
1,2為碰撞前速度，𝜔⃑⃑ 1,2為碰撞前 角速度，𝑛⃑ 為連心線的單位向量，𝑣
𝑐為接觸點相對速度，𝑛⃑和
𝑣 𝑐的夾角為碰撞角度 。
。
圖 3-2
𝜔 ⃑⃑
1𝑣
1𝑣
2𝑛⃑
𝜔 ⃑⃑
2𝑣
𝑐𝛾
顆粒 1
顆粒 2
54
圖 3-3 分子動力學模擬初始情況示意圖:
模擬的初始情形，藍色圓圈是顆粒氣體，紅色圓圈是顆粒布朗物體，邊界為週期性邊
界，一開始顆粒氣體的x、
y
方向速度為常態分布，其顆粒溫度(能量)為T / 2
。初始情況顆粒氣體的角速度為 0，而顆粒布朗物體一開始x、
y
方向速度和角速度都為 0。圖 3-3
顆粒氣體 顆粒布朗物體
週期性邊界
𝑣𝑥=0 𝑣𝑦=0 𝜔 =0
𝐸𝑥=𝑇/2 𝐸𝑦=𝑇/2 𝐸𝜔=0
55
圖 3-4 顆粒氣體和顆粒布朗物體速度分布圖
(a)顆粒氣體(b)顆粒布朗物體之三個自由度速度的分布機率，紅色的點為x方向速
度，黑色的點為
y
方向速度，綠色的點為角速度。圖 3-4
機 率
機 率
速度(公分/秒)，角速度(弧度/秒) (a)
(b)
速度(公分/秒)，角速度(弧度/秒)
𝑽𝒙 𝑽𝒚 𝝎
𝑽𝒙 𝑽𝒚 𝝎
56
圖 3-5 顆粒氣體和顆粒布朗物體能量分布圖:
(a)顆粒氣體(b)顆粒布朗物體之三個自由度能量的分布機率，其中Ex是x方向動
能，Ey是
y
方向動能、E是轉動動能。圖 3-5
機 率
機 率
動能(奈焦耳) (a)
(b)
動能(奈焦耳)
𝑬𝒙 𝑬𝒚 𝑬𝝎
𝑬𝒙 𝑬𝒚 𝑬𝝎
57
58
59
表 3-1 單位對應表
模擬的時間跟真實時間的對應表，模擬總時間和溫度調整的時間間隔以及數據截取 的時間間隔，經過換算得出真實物理時間。
表 3-1
時間對應表 模擬單位 真實物理單位
長度 1 0.1(公分)
單元分區寬度 1.1 0.11(公分) 質量 1 4.47x10
-3(克)
速度 1 0.2(公分/秒)
溫度(能量) 1 1.78x10
-5(微焦耳)
時間 1 0.5(秒)
每ㄧ步時間間隔 10
-55x10
-6(秒)
溫度調整時間間隔 10
-45x10
-5(秒)
數據截取時間間隔 10
-10.05(秒)
模擬總時間 10
95000(秒)
60
表 3-2 模擬參數表和結果。
用實驗量到的參數在分子動力學模擬得到的顆粒氣體和顆粒布朗物體的平均能量。
表 3-2
參數 模擬單位
N
:粒子數量 70a:顆粒布朗物體半徑 8
m:顆粒布朗物體質量 50
T
:溫度 1000(17.8 奈焦耳)
1:顆粒氣體-顆粒布朗物體垂直碰撞係數 0.65
2:顆粒氣體-顆粒氣體垂直碰撞係數 0.65q
:轉動慣量係數 0.5u
:摩擦係數 0.2結果(顆粒氣體)
E
x :x方向平均動能 8.88
0.01 奈焦耳
E
y :y方向平均動能 8.92
0.01 奈焦耳
E
:轉動動能 0.39
0.01 奈焦耳結果(顆粒布朗物體)
E
x :x方向平均動能 8.80
0.13 奈焦耳
E
y :y方向平均動能 8.90
0.13 奈焦耳
E
:轉動動能 0.42
0.01 奈焦耳61
參考資料
[1]R. A. L. Jones, Soft Condensed Matter, Oxford University Press(2002).
[2]Stephen J. Blundell and Katherine M.Blundell, Thermal Physics, Oxford University Press(2010).
[3]賈魯強,黎璧賢,漫談顆粒體物理,物理雙月刊,23卷4期,503 (2001).
[4]T. Wang, K. To, Granular gas in a vibrating box, Chin. J. Phys. 45, 675 (2007).
[5]陳文楠, An experimental study of the motion of a quasi two-dimensional granular gas, 國立清 華大學物理所碩士論文(2008).
[6]W. Chen, K. To, Unusual diffusion in a quasi-two-dimensional granular gas, Phys. Rev. E 80, 061305 (2009).
[7]K. To, Boltzmann distribution in a nonequilibrium steady state: Measuring local potential by
granular Brownian particles, Phys. Rev. E 89, 062111 (2014).
[8]J. S. van Zon, F. C. MacKintosh, Velocity distributions in dilute granular systems, Phys. Rev. E
72, 051301 (2005).
[9]W. Losert, D.G.W. Cooper, J. Delour, A. Kudrolli, J.P. Gollub, Velocity statistics in excited
granular media, Chaos 9, 682 (1999).
[10]F. Rouyer, N. Menon, Velocity fluctuations in a homogeneous 2D granular gas in steady state, Phys. Rev. Lett. 85, 3676 (2000).
[11]J.S. Olafsen, J.S. Urbach, Velocity distributions and density fluctuations in a granular gas, Phys.
Rev. E 60, R2468 (1999).
[12]G.W. Baxter, J.S. Olafsen, The temperature of a vibrated granular gas, Granular Matter 9, 135 (2007).
[13]K. Feitosa, N. Menon, Breakdown of energy equipartition in a 2D binary vibrated granular gas, Phys. Rev. Lett. 88, 198301 (2002).
62
[14]K. Nichol, K.E. Daniels, Equipartition of rotational and translational energy in a dense granular
gas, Phys. Rev. Lett. 108, 018001 (2012).
[15]S. McNamara, S. Luding, Energy nonequipartition in systems of inelastic, rough spheres, Phys.
Rev. E 58, 2247 (1998).
[16]D. C. Rapaport, The Art of Molecular Dynamics Simulation - 2nd ed., Cambridge University Press, (2004).
[17]S. Luding, Granular materials under vibration: Simulations of rotating spheres, Phys. Rev. E
52, 4442 (1995).
[18]A. Barrat, E. Trizac, Molecular dynamics simulations of vibrated granular gases, Phys. Rev. E
66, 051303 (2002).
[19]O. Herbst, R. Cafiero, A. Zippelius, H.J. Herrmann, S. Luding, A driven two-dimensional
granular gas with Coulomb friction, Phys. Fluids 17, 107102 (2005).
[20]A. Barrat, E. Trizac, Lack of energy equipartition in homogeneous heated binary granular
mixtures, Granular Matter 4, 57 (2002).
[21]N.V. Brilliantov, T. Po¨schel, W. T. Kranz, A. Zippelius, Translations and Rotations Are
Correlated in Granular Gases, Phys. Rev. Lett. 98, 128001 (2007).
63
附錄:程式
gif2xyat
#!/bin/bash
# where : gif2xyat.s
# usage : gif2xyat.s a.gif [ARGUMENTS]
# parameter:
# d : the diameter of one ball
# t1 : threshold
# v : verbose
# sw : smooth width
# This program scans the balls for each tif image in the filelist
# require : IDL, find_chain.pro, Eric Weeks's liberary(in the directory "lib") pro=/ramdisk/gif2xyat.$$; run=/ramdisk/idlrun.$$
#pro=/ramdisk/gif2xyat.$$; run=/ramdisk/idlrun.$$
libpath='/home/user/two_side_gif'
mox='0';moy='0';d='12'; t1='120';t2='120';sw2=2;v=";";v2=";"; sw=1; wt=1; br=9;
fn="$1"; shift
cat > $pro << IDLPRO
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
64
; start the main program fn='$fn'
sx=info.dimensions[0] & sy=info.dimensions[1]
for i=0, info.num_images-1 do begin read_gif, fn, im, MULTIPLE=i
$v window, 0, xsize=info.dimensions[0], ysize=info.dimensions[1]
$v tvscl,im
65
cc =smooth(c,$br) idd=where(cc eq 0) im[idd]=0
a=smooth(hist_equal(im,minv=$t2),$sw2)
$v2 window, 0, xsize=info.dimensions[0], ysize=info.dimensions[1]
$v2 tvscl, im
;;;find orientation
rrmax=max(dr,imax,subscript_min=imin) ;& print, imax, dr[imax], imin, dr[imin]
dx=mx1[imax]-mx1[imin] & dy=my1[imax]-my1[imin] & angle=atan(dy,dx) print, format='(%"%f\t%f\t%f\t%d")',mx,my,angle,i
cat > $run << IDLRUN
!PATH=!PATH+":$libpath/lib"
.run $pro
66
exit IDLRUN
idl $run 2> gif2xyat.err rm $pro; rm $run
do.060302.s
#!/bin/bash
# script : do.060302.s
dir=`echo $0 | cut -d. -f3`; [[ -d $dir ]] || mkdir $dir;
cp $0 ~/tools/gif2xyat.s ~/tools/cut.s ~/tools/eng.s $dir; cd $dir;
cat > parameter.txt << PAR type = disk
p2cm = 0.1075 cm/pixel fps = 100 #
tinterval = 0.5 seconds PAR
for ((i=0;i<5000;i++)); do
#for ((i=10;i<11;i++)); do for ((j=0;j<7;j++)); do
t=`echo $i $j | awk '{printf("%05d", $1*8+$2+1)}'`
a=`echo $i $j | awk '{printf("%d", $1*8+$2+1)}'`
[[ $a -gt 5000 ]] && break;
67
function posi { cd $dir;
for f in xyat/0[0-9][0-9][0-9][0-9]; do cat $f | mean.s | cut -f1,2,3;
done > xya.all
gp -f xya.all -wd -x 0 501 -y 0 480 -a1 -ntic -nk -png xya.all.png -s 0.8 1;
cd ../; }
function uvwt { cd $dir;
w=`awk
'BEGIN{pi=4*atan2(1,1)}NR>1{dw=$3-x0;if(dw>4)o-=2*pi;if(dw<-4)o+=2*pi;x0=$3;print NR,$3+o}' $f | linfit.s | cut -d' ' -f2 | awk -v t=$time '{print $1/t}'`;
echo $u $v $w $t done > uvwt
cd ../; }
function engt { cd $dir;
68
cat uvwt | eng.s m 59 I 155 > eng cd ../; }
function grav { cd $dir;
gp -g 3 uvwt uvwt uvwt -u 4:1 -u2 4:2 -u3 4:3 -nk -png t.uvw.png cut -d' ' -f1 uvwt | f7.s nb 30 > u.f7
cut -d' ' -f2 uvwt | f7.s nb 30 > v.f7 cut -d' ' -f3 uvwt | f7.s nb 30 > w.f7 gp -g 3 [uvw].f7 -ylog -png uvw.f7.png cd ../; }
// // purpose : simulate balls bouncing in thin vertically moving box // // usage : m n20 t1e6 T100 M100 Ba 20 20 pf | tee f.rv | Show2d.s
// // compile : g++ `itpp-config --cflags-opt` e.c `itpp-config --libs-opt` -o e // // remarks : f.in needed
//
//////////////////////////////////////////////////////////////////////////////
#include <itpp/itbase.h>
#include <itpp/stat/misc_stat.h>
using namespace itpp;
using namespace std;
#include "m.h"
double max_overlap, max_bondLen, min_z;
69
void SetParameters(int argc, char** argv) { int i;
stepLimit=10; eachSnap=1; ATstep=100; Temp=1; radius=5; RandSeed=0;
waitStep=0; decay=0;
rCut=pow(2.0,1.0/6.0); rCutSq=sqr(rCut);
apha=1;aphat=1;apha2=1;aphat2=1;mass=1;qfac=0.5;mui=10;
rNebrShell=0.3; fwall=0.99; fball=1; dt=1e-4; gravity=0;
Out=1; Ka=1e5; Kb=fball*Ka;
BBB.x=20; BBB.y=20; BBB.z=3; tpTemp=1; bsTemp=1; Temp=1; RTemp=1;
bondLen=1.5;
GetNameList(argc, argv);
for(i=1;i<argc;i++) {
70
RNG.setup(0,Temp);//RNG.setup(0,RTemp);
tpRNG.setup(0,tpTemp); bsRNG.setup(0,bsTemp);
tpExpRNG.setup(tpTemp); bsExpRNG.setup(bsTemp);
Set(Box,BBB.x,BBB.y,BBB.z);
///////////////////////////////////////
theta=asin(2*0.5/(radius-0.5));cn=(int)(2*M_PI/theta)+1;n=cn+tn;
///////////////////////////////////////
nMax=(int)(2*Box(0)*Box(1)); if (n>nMax) n=nMax;
Set(region,Box(0),Box(1),Box(2)+2*vibAmp); Bx2=Box(0)/2; By2=Box(1)/2;
cells=to_ivec(region*1/(rCut+rNebrShell)); cells(2)=1;
nCell=cells(0)*cells(1);
invWid=elem_div(to_vec(cells), region); invWid(2)=1;
nxtBead.set_size(n,false); cellList.set_size(nCell,false);
nebrTabMax=n*20; nebrTab.set_size(2*nebrTabMax,false);
vOff(0)="0 0 0"; vOff(1)="1 0 0"; vOff(2)="1 1 0";
vOff(3)="0 1 0"; vOff(4)="-1 1 0";
strcpy(ParFn,prefix); strcat(ParFn,".par"); ParFile.open(ParFn);
strcpy(XyzFn,prefix); strcat(XyzFn,".Xyz"); XyzFile.open(XyzFn);
strcpy(UvwFn,prefix); strcat(UvwFn,".Uvw"); UvwFile.open(UvwFn);
strcpy(OutFn,prefix); strcat(OutFn,".Out"); OutFile.open(OutFn);
}
void SetupJob()
71
{
GlobalRNG_reset();
if(RandSeed==0)RandSeed=RAND.random_int(); GlobalRNG_reset(RandSeed);
rvBf = new double[3*n]; nBuf=3*n*(sizeof (double));
Init_r(); Init_va(); Init_rsxy(); c.set_size(n);
nebrNow=1; SetBase;
nebrNow=1; SetBase;