感謝鄭伊芳助教翻譯課本習題
問題8.1: Lennard-Jones相互作用的定性性質
寫一個簡短的程序，或者使用一個圖形軟件包來繪製的Lennard-Jones勢（8.1）和相應的力的大小:
u (r)在r = 0.8σ的時候，值是多少？如果r下降到r =0.72σ，減少了10％，u 會增加多少？在r = 2.5σ時，r 的值是多少？在r等於多少的力等於零？
問題8.2: 平衡的方法
1. 在Program md中，被併入DATA語句的初始配置對應到N= 16的粒子，而那些粒子要透過一個線性尺寸L = 6的方形單位與Lennard-Jones勢有相互作用。檢查所有粒子的x和y坐標坐落於0和6之間，設置ΔT= 0.01並運行程式以確保它正常運作。總能量應大約被conserved，而且所有16個粒子的軌跡應在屏幕上被看到。
2. 假設在t= 0時，0≤X≤6的限制被除去，在矩形的單元為Lx=12和Ly=6裡粒子可以自由的移動。將這個變化併入你的程序裡並且觀察粒子的軌跡。系統是否隨著時間的增長變的更隨機或更少隨機？
3. 計算n(t)，在該單元格的左半邊中的粒子的數量，和繪製作為一個t函數的n(t)圖。n(t)的定性行為是什麼？另外計算n（t）的時間平均，並繪製t的函數。左半粒子在系統已達到平衡後，他的平均值是多少？用你在問題7.1發現的結果，跟你的定性結果做比較，如果你做了一個分子N = 64的粒子動力學模擬，達到平衡的方法會得到更好的定義嗎？
問題8.3:
Verlet算法的測試
1. 對應於N= 11粒子以相同的速度往相同的方向移動，考慮以下的初始條件（參見圖8.4）。選擇Lx=Ly=10和Δt= 0.01。
FOR i = 1 to N
LET x(i) = Lx/2
LET y(i) = (i -0.5)*Ly/N
LET vx(i)
= 1
LET vy(i)
= 0
NEXT i
系統最終會達到平衡嗎？為什麼或者為什麼不？
圖8.4: 一個特殊的初始條件的範例;箭頭表示幅值和每個粒子的速度方向。
3. 修改Program md，使該程序在一段預定的時間內運行。使用包含DATA語句的初始條件，並保存所有粒子在t =0.5的位置和速度在一個文件中，然後考慮時間倒轉的過程，例如如果時間的方向被逆轉，會產生的運動。這種倒轉等於讓所有粒子v →−v，粒子會返回到原來的位置嗎？假如把速度反向，一會兒後會發生什麼情況？如果你選擇了一個較小的值Δt，會發生什麼事？
4. 使用包含DATA語句的初始條件，但讓Lx = 12，並保存所有粒子在t = 2.5的位置和速度在一個文件中，然後考慮時間倒轉的過程，粒子會返回到原來的狀態嗎？假如把速度反向，在t = 5而不是在t = 2.5的時候，會發生什麼事？
5. 從你的結果中，你可以得到什麼結論有關軌跡的混亂性質？計算出的軌跡是否和“真”的軌跡相同？
問題8.4: Verlet算法的測試
1. 一個必要檢查分子動力學的程序是總能量與所需的精度是守恆的，使用Program md和確定Δt的值所需的總能量是否有跟在t = 2時間間隔上一個給定的精度互相守恆。一種標準是計算ΔEmax(t)時，最大的差|E(t)−E(0)|，在時間間隔t，其中，E(0)是初始的總能量，而E(t)是時間t的總能量。請確認是否當ΔEmax(t)下降時，Δt變小為了固定t，如果你的程序正常工作，那麼ΔEmax(t)應該會減少為(Δt)2。
2. 由於週期邊界條件的使用，中央單元中的所有點是等效的，該系統是平移不變。解釋為什麼Program md應保護總線性動量，浮點誤差和截斷誤差與有限差分算法可以導致total linear
momentum漂移。編程錯誤也可能因動量守恆的檢查被檢測到。因此，它是一個，在規律的時間內監測total linear
momentum的好方法，也是，假如必要時，重置到零的好方法。SUB check_momentum，列在下面，應該被called在SUB initial獲得初始配置後，也在固定時間間隔的程序的主循環，例如，每100 -
1000個時間步長。Program md在ΔT= 0.01時，對total linear momentum的保存效果為多好？
SUB
check_momentum
DECLARE PUBLIC vx(),vy(),N
LET vxsum = 0
LET vysum = 0
! compute
total center of mass velocity (momentum)
FOR i = 1 to N
LET vxsum = vxsum + vx(i)
LET vysum = vysum + vy(i)
NEXT i
LET vxcm = vxsum/N
LET vycm = vysum/N
FOR i = 1 to N
LET vx(i) = vx(i) -vxcm
LET vy(i) = vy(i) -vycm
NEXT i
END SUB
3. 監控程式保存總能量到多好的程度的另一種方法是分析時間序列E(t)，在一條直線上利用最小二乘擬合E(t)，直線的斜率可以解釋為漂移和root mean square deviation從直線可以被解釋為噪聲（第7.5節中的符號σy）。漂移和噪聲如何依賴Δt，在一個固定的時間間隔t裡？Most research
applications conserve the energy to within 1 part in 103 or 104 or better over
the duration of the run.(看不懂)
4. 考慮高階算法之一，在附錄5A裡曾被討論到的牛頓運動方程式的解決方法，你能使用Verlet的算法去選擇一個較大的Δt值來達到同樣程度的節約能源嗎？總能量是否出現波動和/或到最後漂移？
問題8.5: 定性性質的液體和氣體
1. 對於這個問題，請使用存儲在Program md中DATA語句裡的初始配置（如果你不希望輸入初始條件，請參見問題8.11，來幫助你想出你自己的）。對於這個初始條件N = 16和Lx= LY =6，什麼是較低的密度？什麼是系統的初始能量？選擇Δt = 0.01並且模擬至少500個時間步長或t = 5。用理想氣體的質與你對P的估計做比較。
2. 修改你的程式使溫度和壓力的的瞬時值在系統達到平衡前不會累積。什麼是你對平衡的標準? 其中有一個標準是在有限的時間間隔計算T和P的平均值，並檢查這些平均值不隨時間漂移。
3. 開始仿真的方法之一是保存從較早的運行中使用的位置。最簡單獲得初始條件對應於不同的密度但又具有相同值N的方式是重新調整粒子位置和單元格的線性尺寸的比例。下面的代碼顯示了這樣做的方法之一:
LET rscale = 0.95
FOR i = 1 to N
LET x(i) = rscale*x(i)
LET y(i) = rscale*y(i)
NEXT i
LET Lx = rscale*Lx
LET Ly = rscale*Ly
LET area = Lx*Ly
上面的代碼納入到你的程式在一個單獨的子程序中。你會怎麼期望P和T的改變，當系統被壓縮的時候？使用在（a）部分中相同的初始結構並且確定平均溫度和壓力，在密度ρ = 16/(5.7)2 (平方) ≈ 0.49的時候（重新定標= 0.95）。比較你得到的P和理想氣體的結果。如果你選擇重新調整的質是更小的，你認為會發生什麼事？儲存你的最終配置模擬在一個文件中，並用它作為另一個運行ρ = 16/(0.95 × 5.7) 2 (平方) ≈
0.55的初始條件。決定這種更高的密度下的T和P，並保存最終配置。
問題8.6: 速度和流速分佈
1. 寫一個子程序來計算平衡機率P(v)Δv，其中一個粒子具有v 和v+Δv之間的速度。要這樣做，預估最大速度的值vmax，你需要放入箱中，選擇箱的寬度dv = vmax/nbin的，其中nbin是箱的總數。下面的代碼說明了一種把速度到其應有的箱(bin):
LET v2 = vx(i)*vx(i) + vy(i)*vy(i)
LET v = sqr(v2)
LET ibin = truncate(v/dv,0) + 1
IF ibin > nbin then LET ibin = nbin
LET prob(ibin) = prob(ibin) + 1
數組元素的prob(ibin)記錄粒子的速度對應於ibin的次數。雖然這是低效的，在每一個時間步長至少為100個時間步長後決定prob(ibin)。規格化prob(ibin)除以通過的粒子的數量和時間步長的數量。使用初始配置存儲在DATA語句在Program md對於這個問題。選擇nbin=50。
2. 畫出概率密度P(v)相對於v，P(v) 的定性形式是什麼？v最有可能的值是多少？ P(v)的大致寬度是多少？和理論形式比較你的測量結果（在兩個方面）
速度分佈的形式（8.11）被稱為Maxwell-Boltzmann分佈。
3. 決定速度x和y分量的概率密度，確保你有區分正和負價值。x和y速度分量最可能的值是多少？其平均值是多少？繪製概率密度P(vx) 相對於 vx 和 P(vy) 相對於vy。透過繪製平均值[P(vx = w) + P(vy
= w)]/2 相對於 w，你可以發現更好的結果。P(v)的定性形式是什麼？它和速度的概率密度相同嗎？
問題8.7: 溫度和壓力的能量的依賴關係
1. 我們已經在問題8.5看到，總能量由初始條件決定，且溫度是由只有在該系統已達到熱平衡才能導出的量。因為這樣，在一個特定的溫度下要研究系統是很難的，藉由重新調整系統的速度，平均溫度可以被改為所需的溫度，但我們必須小心，不要一下太快增加速度。為什麼? 從較早的模擬Lx = Ly = 6 和 N
= 16中，選擇初始配置的系統為一個平衡的配置，並決定T(E)平均溫度的能量依賴性，在範圍T = 1.0到T = 1.2。由所需的因子重新調整速度，平衡溫度是否提高到所需的值？如果沒有，你需要再次重新調整速度。在一般情況下，所要求的溫度是由一系列重新調整速度後才能達成，並且要在夠長時間裡，如此一來才能在調整期間讓系統維持接近平衡。
2. 使用你在(a)部分發現的數據T(E)繪製總能量E的函數，T是單調遞增E的函數嗎？從potential energy去預估對CV的貢獻，由於potential energy的關係，有多少百分比的貢獻是給熱容量？為什麼精確測定的CV難以實現？
3. 在我們的分子動力學的模擬裡，總能量是固定的，但可能會發生波動的動能和勢能。另一種確定CV的方法是它與波動的動能有相關性，（在？？章，我們找到CV與總能量的波動有關，在constant T, V,N ensemble）。可以被這麼表示（參見Ray和Graben）
我們要怎樣才能為在高密度的粒子系統產生一個典型的初始條件？如果我們單純將粒子隨機放在中央cell，會發生什麼情況？現在的問題是，如果系統是密集的，一些粒子會彼此非常接近，並對彼此產生一個非常大的排斥力F。由於條件(F/m)(Δt)2 << σ必須被適用的有限差分方法(finite difference method)來滿足，粒子的隨機放置是不實用的。但是，如果一個虛構的曳力速度與速度的平方成比例被引入到平衡系統，要使用這樣的初始條件是有可能的。這種力量的影響是那些粒子速度的減幅，而那些粒子會變得過大，因為有很大的力施加在他們的身上。
問題8.8:
圖8.5: 在一個三角晶格每個粒子有6個最近的鄰居。
從圖8.5可以看出三角晶格的性質，每個粒子有6個最近的鄰居。雖然它是可以選擇三角晶格的中央cell為一個菱形，選擇單元為矩形的會是更方便，我們採取cell的線性尺寸分別為Lx 和 Ly =√3Lx/2。為了簡單起見，我們假設√N是一個整數，因此，在水平和垂直方向上的晶格間距各別為ax = Lx/√N 和 ay = Ly/√N，在每行中的晶格位置被前一行取代1/2ax。下面的代碼生成三角晶格:
LET i
= 0
FOR col = 1 to nx
FOR row = 1 to ny
LET i = i + 1
IF (mod(row,2) = 0) then
LET x(i) = (col
-0.75)*ax
ELSE
LET x(i) = (col
-0.25)*ax
END IF
LET y(i) = (row -0.5)*ay
NEXT row
NEXT col
寫一個程式來計算透過Lennard-Jones勢相互作用在N粒子的一個系統裡每個粒子的勢能，考慮三角形和正方形晶格兩者，選擇正方晶格的線性尺寸為L = _LxLy，使兩個晶格具有相同的密度。選擇N =36，並確定能量為Lx = 5和 Lx = 7，每一種情況下系統的密度為多少？ 你E/N的結果是取決於晶格的大小嗎？哪種晶格對稱具有較低的能量？用三角晶格把粒子包的更加接近的能力來解釋你的結果。
問題8.9: 固態和熔融
1. 選擇N = 16，為Lx = 4，和Ly =√3Lx/2，並在粒子上放置三角晶格。給每個粒子的初始速度為零，系統的總能量是多少？做一個模擬，並以時間的函數測量溫度和壓力，該系統是否仍是固體？
3. 從（b）部分選擇初始配置為一個平衡配置，但藉由兩個因素增加其動能，新的總能量是多少？描述粒子運動的定性行為，該系統的平衡溫度和壓力是多少？達到平衡後，以相同的方式重新調整速度再次增加溫度，重複此重新縮放並測量P(T)和E(T)的幾種不同的溫度。
4. 使用你在(c)部分的結果，並設計函數為T的E(T) − E(0)和P(T)，E(T) − E(0)的差異與T有成比例嗎?什麼是和諧固體的平均勢能？其熱容量是多少？
5. 從（b）部分選擇一個平衡配置，並藉由重新縮放Lx, Ly和1.1因素去做協調的粒子減少密度。軌跡的本質是什麼？
減小系統的密度，直到系統熔化。你的熔化質量標準是什麼？
問題8.10: 亞穩定性
1.
我們可以建造一個亞穩態，藉由把粒子放在一個正方形的cell裡，其形狀與最低能量狀態相對應的三角晶格不符。修改SUB initial使粒子的位置形成一個方形的晶格，如果初始速度設置為零，當你運行程式的時候會發生什麼事？選擇N = 64和Lx = Ly = 9。
3. 如果時間允許，請重複(b)部分與vmax = 0.1