2.1.2 麥克斯韋-波爾茲曼能量分佈
目的：確認系統中原子的動能分佈與給定溫度下的麥克斯韋分佈相對應。
對於本練習，您需要一個文件，其中包含模擬一瞬間每個原子的速度。您將使用泛素蛋白平衡的最後一幀的相對應文件。生成此文件的模擬與您在目錄1-3-box /中執行的模擬完全相同，但此處使用的時間步長為1fs。您可以將用於該模擬的配置文件與該目錄的example-output /中的配置文件進行比較。正如第一單元中討論的那樣，在任何生產模擬中，都應將rigidBonds用於水，因為水分子已經參數化為剛性分子。但是，對於麥克斯韋-波爾茲曼能量分佈的說明，力場的影響可以忽略不計。為簡單起見，我們沒有使用選項rigidBonds。在這裡，我們為您提供文件ubq_wb_eq_1fs.restart.vel。或者，可以使用任何速度文件，只要模擬中未使用約束。
1、現在將加載結構文件。首先，打開Molecule File Browser（分子文件瀏覽器），File→ New Molecule菜單項打開。瀏覽並加載公共目錄中的文件ubq_wb.psf。
2、現在，窗口頂部的菜單應顯示為3：ubq_wb.psf。 （ps：教程裡默認你的VMD一直沒有關過，這在學習過程是不可能的，因為關掉重新開始還是0，而且，0,1,2,3,4等等沒有關係）這樣可以確保將您加載的下一個文件添加到該分子（分子ID 3）。現在，瀏覽目錄2-2-maxwell /中的文件ubq_wb_eq_1fs.restart.vel。這次您需要選擇文件類型。在“確定文件類型：”下拉菜單中，選擇“ NAMD Binary Coordinate”，然後再次單擊“Load”。關閉分子文件瀏覽器。
您加載的分子看起來很可怕！那是因為VMD正在讀取速度，就好像它們是坐標一樣。它的外觀無關緊要，因為我們想要的是使VMD編寫一個有用的文件：該文件包含每個原子的質量和速度。
3、在VMD TkCon窗口中，通過鍵入cd ../2-2-maxwell/導航到2-2-maxwell /目錄。
4、首先，您將創建一個包含系統中所有原子的選擇。通過在TkCon窗口中鍵入以下內容來進行設置：
set all [atomselect top all]
5、現在，您將打開文件energy.dat進行寫入：
set fil [open energy.dat w]
6、對於系統中的每個原子，計算動能，並將其寫入文件
foreach m [$ all get mass] v [$ all get {x y z}] {puts $ fil [expr 0.5 * $ m * [vecdot $v $v]]}
7、關閉文件：close $ fil
8、在VMD外部，使用寫字板查看新的energy.dat文件。打開文件時，請確保下拉菜單“文件類型”設置為“所有文件”。它看起來應該像下面的列表（儘管值不一定相同）：
1.26764386127
0.189726622306
0.268275447408
...
9、退出寫字板。
10、如果您選擇不執行上述步驟，則有一個VMD腳本可以復制它們並創建文件energy.dat。在“ TkCon”窗口中，鍵入：
source get_energy.tcl
這將獲取腳本get_energy.tcl，並生成energy.dat文件。
現在，您有一個原始數據文件，可用於將獲得的能量分佈擬合到Maxwell-Boltzmann分佈。您可以從該擬合中獲得分佈的溫度。在本節的其餘部分，我們將展示如何使用Excel進行此操作。如果您更熟悉其他工具，則可以改用它。
11、退出VMD。
12、在Excel中，單擊文件→ 打開...，導航到目錄2-2-maxwell /，然後選擇文件energy.dat。確保下拉菜單文件類型設置為所有文件。
13、出現“文本導入嚮導”時，為文件類型選擇“分隔符”，然後單擊“完成”。能量數據將出現在列A中。
我們將構建能量值的直方圖，並對其擬合MaxwellBoltzmann分佈，以找到模擬的平均溫度。
14、構造的直方圖一列柱標籤。在單元格B1中，輸入數字0，在單元格B2中，輸入公式“= B1 + 0.1”，然後按Enter鍵。單擊單元格B2，單擊並按住其右下角，然後向下拖動至第71行並釋放。您應該具有一個從0到7填充B列的值列表。這些是我們的直方圖柱。
15、現在，我們將使用數據分析工具來創建直方圖數據。點擊工具→ 數據分析...。如果該選項不存在，則需要通過單擊工具→ 添加插件...並從列表中選擇分析工具庫來添加數據分析工具。
（ps：這段我不是這樣操作的，因為Excel沒有“工具”這一項）
（ps：假如沒有數據分析這一項，應該按照：文件——選項——加載項——分析工具庫）
16、從數據分析窗口中選擇直方圖，然後單擊確定。輸入以下單元格範圍，然後單擊“確定”：
輸入範圍：$ A $ 1：$ A $ 7051
Bin範圍：$ B $ 1：$ B $ 71
將使用您的能量柱以及給定能量柱中發生的能量測量數量創建一個新表。我們可以繪製這些值，但是我們寧願繪製歸一化頻率。
要獲得歸一化的頻率，我們將獲得直方圖下方的面積，然後將每個頻率除以該面積。
17、在新工作表的單元格C2中，鍵入公式”= B2 * 0.1”。然後單擊並按住單元格C2的右下角，然後向下拖動到第72行以填寫每個bin的面積值。
18、在單元格D2中，輸入公式”= sum（C2：C72）”。這是直方圖下的區域。
19、通過在單元格E2中鍵入以下公式來創建歸一化頻率：”= B2 / D $ 2”。點擊enter。像之前對C列所做的那樣，將此輸入向下擴展到第72行。
20、通過單擊插入→ 圖表...，選擇柱狀圖作為圖表類型，選擇左上角的第一個圖表子類型，然後單擊下一步，來繪製數據的直方圖。
21、在步驟2中，單擊“系列”選項卡並刪除所有現有系列。然後通過單擊添加按鈕添加您自己的。單擊“值：”字段中的紅色按鈕，選擇E列（行2-72）中的每個數據值，然後單擊紅色按鈕完成操作。對類別（X）軸標籤執行相同的操作：使用A列。單擊完成。
22、隨時修改圖表。我們建議右鍵單擊直方圖欄，選擇“格式數據系列...”，單擊“選項”選項卡，並將“間隙寬度：”設置為0。
現在，您將使Maxwell-Boltzmann分佈適合此圖，並找出該分佈對應的溫度。 （希望，您要在仿真中達到的溫度！）
23、在單元格F2中，輸入公式：
“=（2 / SQRT（pi（）* G $ 2∧3 ））* SQRT（A2）* EXP（-A2 / G $ 2）”
輸完之後，在單元格G2中輸入值0.5，然後將F2中的公式向下擴展到單元格F72。這將創建Maxwell-Boltzmann分佈，我們將其擬合原始數據。單元格G2中的參數對應於以kcal / mol為單位的k B T。 （這些是NAMD使用的能量單位。）
24、在柱狀圖，在柱狀圖單擊鼠標右鍵，並點擊“選擇數據”。單擊“添加（A）”，生成一個序列名稱為2，系列值為F列第2-72行，單擊確定。 （ps：跟原英文的效果是一樣，我按照中文版本EXCEL做的圖，如下圖會出現兩個顏色的柱狀圖）
25、右鍵單擊剛創建的紅色直方圖條，選擇“圖表類型...”，選擇“線”作為圖表類型，然後選擇左上角的第一個“圖表”子類型。單擊確定。
您可以通過右鍵單擊直方圖中的線並選擇“設置數據序列格式...”來更改Maxwell-Boltzmann擬合的線寬。
我們繪製了Maxwell-Boltzmann表達式，但是擬合不是最佳的。確實，如果我們更改單元格G2中k B T的值，我們將看到擬合度變化。現在，我們將確定產生最佳擬合的k B T的值。
26、在單元格H2，鍵入公式“=（F2-E2）∧ 2”.然後單擊並按住細胞H2的lowerright角和向下拖動到行72填寫為每個段的誤差值。
27、在單元格I2中鍵入= sum（H2：H72）。這是我們擬合的總平方誤差，也是我們尋求最小化的誤差。我們將使用Excel中“規劃求解器”。
28、個單擊工具→ 求解...。如果該選項不存在，則需要通過單擊“工具” →“ 添加項...”並從列表中選擇“求解器加載項” 來添加“求解器”工具。
（ps：教程中的excel和我版本不一樣，我的是點擊“數據”——“規劃求解”）
29、在“規劃求解參數”窗口中，設置以下單元格範圍，然後單擊解決：
設置目標單元格：$ I $ 2
等於：最小值
通過更改可變單元格：$ G $ 2
規劃求解器將找到k B T（單元格G2）的值，該值最合適地將Maxwell-Boltzmann分佈與我們的原始能量數據進行擬合。你的圖應該像圖10。獲取溫度T為這種分佈。您的結果有望接近310 K！
您的值與預期價值之間的差異是由於採樣不佳以及蛋白質的有限大小所致；後者意味著即使在完美採樣的情況下，也只能以大約±10 K的精度“測量”溫度。
30、如果需要，保存在Excel中創建的文件，然後關閉Excel。
沒有留言:
張貼留言