2.2.1 熱擴散
目標:模擬泛素分子的冷卻,並根據模擬結果確定係統的熱擴散率。
在本練習中,您將使用水球中泛素平衡的最後一步,文件ubq_ws_eq.restart.coor。 (ps:用到之前模擬的最後一步的坐標,速度等來重新開始另一個模擬) ,而水泡的其餘部分設置為300K。這樣,您可以確定通過監視系統溫度的變化並將其與理論表達式進行比較來確定熱擴散率。
熱擴散。局部溫度T(~r;T)的時間依賴性由熱擴散方程控制
根據初始和邊界條件,對於半徑為R的球體
其中D是熱擴散係數,Tsim是球體的初始溫度,Tbath是球體邊界的溫度。這種系統的平均溫度取決於
1、在“CMD”窗口中轉到該問題的目錄:
cd ../2-6-heat-diff/
為了使用NAMD的溫度耦合功能,您需要創建一個文件來標記受溫度耦合影響的原子。以下幾行:
tCouple on
tCoupleTemp 200
tCoupleFile ubq_shell.pdb tCoupleCol B
在提供的配置文件ubq_cooling.conf中將啟用此功能,並將耦合到熱浴的原子B列的值設置1 。在下面創建的ubq_shell.pdb的B列中為1.00。
2、在VMD TkCon窗口中,導航到2-6-heat-diff /directory。 cd ../2-6-heat-diff/
3、通過在VMD TkCon窗口中鍵入以下內容,將系統加載到VMD:
mol../common/ubq_ws.psf
mol addfile ../1-2-sphere/ubq_ws_eq.restart.coor
4、選擇系統中的所有原子:
set selALL [atomselect top all]
5、查找系統中心:
set center[measure center $selALL weight mass]
6、找到系統中心的x,y和z坐標,並將它們的值放在變量xmass,ymass和zmass中:
foreach {xmass ymass zmass} $center {break}
7、在外層中選擇原子:
設置shellSel [atomselect top“not(sqr(x- $ xmass)+ sqr(y- $ ymass)+ sqr(z- $ zmass)<= sqr(22))“]
8、將此選擇中的原子的beta參數設置為1.00:
$shellSel set beta 1.00
9、創建一個pdb文件,該文件在beta列中用“ 1.00”標記外層中的原子:
$ selALL writepdb ubq_shell.pdb
10、運行分子動力學模擬
namd2 ubq_cooling.conf > ubq_cooling.log &
11、如果您自己不運行模擬,則需要從提供的日誌文件中獲取溫度數據。從目錄example-output /複製此日誌文件。
copy example-output \ ubq_cooling.log。
如果您自己運行模擬,則日誌文件將由NAMD生成,並且您無需複制它。
12、使用namdstats.tcl來了解系統的溫度如何隨時間變化。在“ VMD TkCon”窗口中,鍵入:
source../2-3-energys/namdstats.tcl
data_time TEMP ubq_cooling.log
TEMP.dat文件將隨著時間的推移寫入系統溫度。
13、現在,您可以退出VMD。
14、在Excel中,單擊文件→ 打開...,導航到目錄2-6-heat-diff /,然後選擇文件TEMP.dat。確保下拉菜單文件類型設置為所有文件。
15、出現“文本導入嚮導”時,為文件類型選擇“分隔符”,然後單擊“下一步”。在“分隔符”類別下,單擊“選項卡”。 “數據預覽”窗口應將您的時間步長和溫度數據分為兩列。單擊完成。時間步長和溫度值將分別出現在列A和B中。
16、在繪製數據之前,我們將修改時間單位,使其以fs為單位,而不是時間步長。在單元格C1中,鍵入公式= A1 * 2,單擊並按住C1的右下角,然後向下拖動到單元格C201。這些值是以fs為單位的仿真時間。
17、通過單擊插入→ 圖表...,選擇XY(散點)作為圖表類型,選擇右下角的最後一個圖表子類型,然後單擊下一步,繪製溫度數據。
18、在步驟2中,單擊“系列”選項卡,刪除任何現有的系列,然後單擊“添加”按鈕。單擊“ X值:”字段中的紅色按鈕,選擇列C中的每個數據值(行1-201),然後單擊紅色按鈕完成操作。對Y值執行相同的操作:使用B列中的數據。單擊“完成”。
19、雙擊y軸並將“最小值:”值更改為200。單擊“確定”。如果需要,可以進一步修改圖表。
20、現在,您將使模擬的溫度擬合理論表達式。
21、在單元格D1中,鍵入或複制並粘貼公式(確保它在一行上):
= 200
+ 66.87 *(EXP(-0.0146 * E $ 1 * C1)+ 0.25 * EXP(-0.25 * 0.0146 * E $ 1 * C1)+ 1/9 * EXP(-1 / 9 * 0.0146 * E $ 1 * C1 )+ 1/16 * EXP(-1 / 16 * 0.0146 * E $ 1 * C1)
+ 1/25 * EXP(-1 / 25 * 0.0146 * E $ 1 * C1)+ 1/36 * EXP(-1 / 36 * 0.0146 * E $ 1 * C1)
+ 1/49 * EXP(-1 / 49 * 0.0146 * E $ 1 * C1)+ 1/64 * EXP(-1 / 64 * 0.0146 * E $ 1 * C1)
+ 1/81 * EXP(-1 / 81 * 0.0146 * E $ 1 * C1)+ 1/100 * EXP(-1 / 100 * 0.0146 * E $ 1 * C1))
單元格E1中的值是熱擴散率,我們將對其進行更改以獲得更好的擬合度。在單元格E1中輸入0.03的初始猜測。然後將D1中的公式向下擴展到單元格D201。
22在圖表中,在溫度數據單擊鼠標右鍵,右擊並點擊“選擇數據”。單擊系列選項卡,並添加具有D列第1-201行的Y值的系列, X值:還是C列。單擊確定。
您可以通過右鍵單擊直方圖中的線並選擇“設置數據序列格式...”來更改適合的線寬。
我們繪製了理論溫度表達式,但是擬合不是最佳的。如果我們更改單元格E1中的熱擴散率,則會看到擬合度變化。現在,我們將確定最佳參數。
23、在小區F1,鍵入公式=(D1-B1)∧ 2.然後單擊並按住小區F1和拖動的lowerright角落向下到行201填補對每個段的誤差值。
24、在單元格G1中鍵入= sum(F1:F201)。這是擬合我們的總平方誤差,也是我們尋求使用規劃求解工具使之最小化的原因。
25、點擊數據→ 規劃求解。如果該選項不存在,則需要通過單擊“工具” →“ 添加項...”並從列表中選擇“求解器加載項” 來添加“求解器”工具。 (ps:文件-選項-加載項-找到規劃求解點擊-轉到-確定)
26、在“規劃求解參數”窗口中,設置以下單元格範圍,然後單擊
解決:
設置目標單元格:$G$1
等於:最小值
通過更改單元格:$E$1
Solver將找到最適合我們原始溫度數據的溫度表達式的熱擴散率。
27、將單元格E1中的參數值乘以0.1得到以cm 2 s -1為單位的熱擴散率。您應該獲得大約0.45 ×10 -2 cm 2 s -1的值您的結果。與水的熱擴散率D =1 .4×10 -3 cm 2 s -1相比如何?
28、計算該系統的導熱性,ķ,由下式:ķ = ρC V d,其中Ç V是比熱(在章節2.1.5計算),d是熱擴散率(假設系統密度ρ為1 g / mL)。
29、如果需要,保存在Excel中創建的文件,然後關閉Excel。
沒有留言:
張貼留言