所謂分子動力學（MD）模擬，是在預先確定條件下（比如溫度，壓力，應力，外力等等）模擬原子和分子運動的一種方法。分子動力學模擬可以用來研究納米尺度下的動力學過程，還可以用來計算相圖、擴散係數和各種響應函數等諸多性質。
QuantumATK可以使分子動力學模擬變得非常簡單：只需要向構型（configuration）添加所需的計算器（calculator）和分子動力學模塊（MolecularDynamics block），就可以開始計算了。對於不同類型的模擬，什麼樣的參數和設置是最合適的，在這裡給出了一些指導方針。本文簡要介紹 QuantumATK 中分子動力學（MD）功能，並且一步步的解釋了如何正確地設置模擬過程以得到想要的結果。
分子動力學模擬本質上是對給定初始構型中原子運動的牛頓方程求數值解。這通常是通過對有限時間步長進行積分計算得到的。原子之間的相互作用（也就是原子間力）可以通過不同的方法（可以是密度泛函理論（DFT）或經典力場）來進行計算。這些力決定了原子的加速度並且使得原子的位置和速度傳遞到下一個時間步長。多次重複這個過程會產生一系列的構型快照，這些快照描述了系統在相空間中的運動軌跡，從而可以進一步從中分析提取出想要的性質。
在設置模擬之前你需要選擇你所感興趣的計算類型。要考慮得因素有：總能量是否守恆（孤立體系中的情況）？模擬體系與熱浴（heat bath）之間的耦合時是否要保持溫度不變？系統是否受到外加壓力或者應力？基於這些考慮我們應該選擇一套合適的模擬參數：時間步長大小， 積分步長的數量（模擬持續的時間），積分算法，初始溫度，約束條件等等。
類似於真實的實驗，進行分子動力學模擬需要一些實證上的知識和經驗。在本例中，你將會了解分子動力學模擬的基本組成部分，以及如何使用 QuantumATK 的分子動力學模塊來進行分子動力學模擬。
將這些技巧應用於一個簡單、著名的測試體系——塊體矽，你將逐步學習到如何選擇時間步長大小，如何控制溫度和壓力，以及如何通過約束條件來固定某些原子。
正如介紹中所述，分子動力學模擬基於的是求解牛頓運動方程。微正則系綜（NVE）是分子動力學模擬中最原始和純粹的形式，其中原子數（N），體積（V）和總能（E）守恆。這些條件對應了一個完全孤立體系，而產生的系綜也被稱為微正則系綜。
在這個章節，你將了解如何運行微正則系綜（NVE）模擬，以及如何通過監測總能（此處應為定值）來評定系統的特性。你將會學習到如何為給定系統選擇合適的時間步長大小。
打開 QuantumATK，創建一個新的項目並命名，打開此項目。點擊工具欄的按鈕啟動 Builder。在 Builder 中，點擊Add ‣ From Database。搜索 Silicon 並添加到 stash 欄中。
初始的矽的單胞的晶格矢量是非正交的。儘管在這種晶胞上可以進行分子動力學模擬，但是使用正交晶胞往往更方便一些（尤其是在晶胞尺寸在模擬過程中發生變化的情況下）。你可以使用 Builder 面板中 Bulk Tools ‣ Supercell。點擊 Conventional 和 Transform 使初始單胞變為常規的立方晶胞。
此外，通常有必要增大矽結構的尺寸以包含更多的原子，因為這樣可以改善所要計算的可觀測量的平均值，並且減少小的模擬晶胞與其周期重複構象相互作用所導致的有限尺寸效應（晶胞矢量的長度最好大於兩倍的勢的相互作用距離）。通過點擊 Bulk Tools ‣ Repeat 選擇 4 x 4 x 4 的重複來增加結構的尺寸。這將得到總共 512 個原子。
接著使用 Builder 窗口右下角的箭頭圖標將結構送到 Script Generator。
接下來，你將設置運行分子動力學計算中將要用到的計算腳本。對於 Si-Si 相互作用，我們將使用 Tersoff 型的勢。實際的分子動力學模擬也將在這裡設置。
在 Script Generator 中
添加一個 New Calculator。
添加 Optimize > MolecularDynamics 模塊。
雙擊 New Calculator 選擇 ATK-Classical 計算器和 Tersoff_Si_1988 勢1)。由於我們不需要保存或者列印初始構型的附加信息，所以不選 Save和 Print。
現在，雙擊 MolecularDynamics 模塊。首先，看一下 Type 下拉列表。這一設置非常重要，它可以讓你選擇使系統處於什麼樣的積分算法和物理系綜下模擬。對於第一個模擬，我們設置為 NVE Velocity Verlet 類型。這個設置選擇了恆定粒子數（N），恆定體積（V）和恆定總能（E）的系綜，並且使用了 velocity verlet 積分算法2)。
設置 Steps 的數量為 1000，Log Interval 為 5。後者的值決定了當前構型的快照多久被寫入硬碟一次。過小得log interval 值會降低模擬的速度，而且會在較大體系計算中產生很大的輸出文件。
在多數分子動力學應用中沒有必要以很高的頻率記錄快照，除非你確實想要分析高頻振蕩。如果你覺得保存快照到硬碟過度的減慢了模擬速度，你可以點掉 Save Trajecotry 複選框來避免模擬過程中快照的寫入硬碟。你可以使用底部的 Save 選項僅保存模擬結束時的軌跡。由於技術原因，這將在保存相同數量的信息的基礎上，大大增快模擬速度。不足的是，當模擬中斷時，你將不能接著中斷後的快照來重啟計算。
鍵入一個合適的軌跡文件名，比如 traj-md-nve-1fs.hdf5。設置 Initial Velocity 為 Maxwell-Boltzmann，選擇 Temperature 為 600 K。這將給每個原子一個隨機的初始速度，對應於 600K 的麥克斯韋-玻爾茲曼分布。另外，選擇相應的複選框，質心的動量將會在分布中被去除從而防止模擬過程中整個系統的漂移。
Time step 在分子動力學模擬中是一個關鍵參數，它決定了數值積分方案的精準度和效率。後續我們會研究它的作用，但是在最開始，我們保持他的默認值 1 fs。
之後，由於軌跡（trajectory）已經被保存，你可以點掉 Save 和 Print 這兩個選項。最終，你的分子動力學設置應該如下圖所示。
點擊 OK 並保存腳本為 md-nve-1fs.py。然後發送腳本到 Job Manager (再次使用箭頭按鈕)並開始模擬。幾分鐘後即可完成模擬。
使用Movie Tool插件，你可以研究系統的溫度、動能、勢能以及總能隨模擬時間變化的函數。同時，原子運動的動畫也將同時展示。選擇LabFloor中的 traj-md-nve-1fs.hdf5 然後點擊 Movie Tool 插件。
如果你查看一下溫度曲線，你就會發現溫度如何迅速下降然後在 300 K 左右進行振蕩。這是因為顯示的溫度是系統的即時溫度，是通過原子的動能以及公式 Ekin(t)=3/2NkBT(t) 計算得出的。由於原子不是孤立的，而是與晶格裡的其他鄰近原子相互作用，他們的初始動能將會部分地轉化成勢能。由於初始構型很接近勢能最小值，接近一半的動能被轉化，導致溫度下降大約50%。實際上，如果你仔細觀察勢能曲線，你會看到相反的情形（勢能迅速升高，然後在一個固定值左右進行振蕩）。
由於你模擬了一個微正則（NVE）系綜，一個重要的監測量是總能的守恆。使用能量曲線圖下面工具欄裡的縮放工具，放大總能曲線。你會看到一個類似於（但可能不會完全一樣）下圖的曲線。
總體上，系統的能量保持在一個相對於平均值小於 10−5量級上的小漲落內。重要的是，平均值是固定的，沒有發現漂移。這些發現表明所選分子動力學（MD）的設置較為恰當。
我們看看如果選擇了一個很大的時間步長大小會發生什麼。修改初始溫度為 5000 K 和時間步長大小為 5 fs，重複上述模擬。你可以再一次使用MolecularDynamics 模塊，也可以將 md-nve-1fs.py 直接拖拽到 Editor 按鈕。記得改一下軌跡（trajectory）文件名。保存這個新修改的腳本為 md-nve-5000K-5fs.py 並運行。
如果你在 Movie Tool 中可視化軌跡模擬的結果，你將發現原子的運動增強。結構甚至會部分融化。如果你再次放大總能曲線，你會觀察到其值顯著的漂移。由於在此溫度下時間步長被選的太大，能量不再變得那麼恆定。
核心要點是，大的時間步長看上去會提高模擬效率，但實際上會增加運動方程數值積分方案的誤差。原子力在一個積分步長內近似守恆的基本假設失效。本質上，所選時間步長應該足夠小以解析原子的最高振動頻率（也就是說，其應該遠小於最小的振動周期）。所以如果你有輕原子（比如氫），那麼通常你需要選擇比只有重原子(比如金)時更小的時間步長。如果在你的計算中有不同的元素，且溫度較高或者原子遠離其平衡構型（即粒子受到很大的力）時，也需要選擇小的時間步長。如果你不知道使用多大的時間步長，對於多數系統你可以選擇 1 fs 作為一個比較安全的開始。通過監測感興趣的條件下的 NVE 模擬總能是否守恆，你可以對大一些的時間步長值進行評估。
另一個需要知道的分子動力學模擬的重要方面是，一個系統需要一定的時間來與所設置的外界溫度和壓力達到平衡。在本例中，系統溫度和總能總的漲落的量級達到穩態用了大約 400 fs。這點很重要，因為任何可觀察量的測量需要在系統達到平衡之後才可以實現（除非你要研究非平衡現象）。在本例中你需要等待至少 400 fs。總的來說，分子動力學模擬的平衡時間最大可以直到幾百納秒（比如在生物分子或者聚合物系統），並且依賴於系統中出現的最長的弛豫時間。所以你應該仔細檢查，所感興趣的可觀察量是否在你的模擬中達到穩態。
在接下來的部分你將學到如何在 NVT 系綜下進行分子動力學模擬，也就是溫度可控的模擬。NVT 系綜也叫作正則系綜，它與微正則系綜 NVE 主要的差別為：系統不再是孤立的，而是可以和一個周圍虛擬的熱浴進行能量交換。
你將使用和上面部分同樣的系統——塊體矽，使用 ATK-Classical 計算器（calculator）以及 Tersoff-potential。
你將以室溫下結構的平衡態為開始。用上述同樣的方法和同樣的模塊（NewCalculator 和 Optimize ‣ MolecularDynamics) 來準備腳本。實際上，如果你先前的腳本還存儲在 ScriptGenerator 裡，你可以重新使用它。
雙擊 MolecularDynamics 模塊打開設置。由於我們要運行 NVT 系綜下的模擬，我們這回需要選擇一個不同的 Type。對於 NVT 模擬，QuantumATK 提供了四種選項：
NVT Berendsen 3)
NVT Nose Hoover4)
Langevin5)
NVT Nose Hoover Chain6)
在 MolecularDynamics 界面工具欄，選擇 NVT Nose Hoover chain 類型。設置步數為 5000 和 log interval 為 20。保存軌跡（trajectory）到文件 traj-nvt-nosehooverchain-300K.hdf5 選擇 Maxwell-Boltzmann 初始速度為 300 K。時間步長默認 1 fs，Reservoir temperature （虛擬熱浴溫度）選擇為 300 K，Thermostat timescale 選擇 100 fs。後面的值決定了系統溫度多快達到虛擬熱浴溫度。我們將會在本部分後面仔細研究這個參數的影響。Final temperature值總是和虛擬熱浴溫度同步，除非你在這裡明確的指定一個不同的值。這個參數可被用來在模擬過程中線性增大或線性減小溫度（參見參考手冊 NVTNoseHooverChain）。Thermostat chain length 這個參數決定了調用多少個並發的調溫器（thermostat）來控制溫度。默認值三可以在多數情況下夠用。只有當你面臨一個強的持久的溫度振蕩時，可以增加該值以獲得更穩定的模擬。
點擊 OK 關閉分子動力學設置，發送腳本到 JobManager 並運行。
模擬結束，軌跡文件將會出現在 Lab floor。使用 Movie Tool 打開。
觀察一下溫度曲線。你會發現在一個短暫的初始平衡時期之後，曲線在所選的模擬熱浴溫度 300 K 附近輕微漲落。實際上，這與之前部分 NVE 模擬的結果看上去沒有太大的不同。與 NVE 模擬不同的是，系統最終平衡溫度不依賴於通過 Maxwell-Boltzmann 分布設定的初始溫度。現在讓我們更細緻的觀察。
打開上一個模擬的 Script Generator。雙擊 MolecularDynamics 模塊改變虛擬熱浴（Reservoir）溫度為 1000 K。Maxwell-Boltzmann 溫度為默認 300 K。最後，將文件名改為 traj-nvt-nosehooverchain-1000K.hdf5。關閉設置，運行腳本。
結果應該與下圖類似。
溫度曲線起始於初始值 300 K（通過 Maxwell-Boltzmann 分布選取）。接著，由於系統與熱浴接觸，溫度迅速提升到虛擬熱浴（reservoir）溫度 1000 K。
你也會注意到總能的增加。與 NVE 模擬相比，這並不表示參數選取不合適，而是反映了系統在與虛擬熱浴交換能量這一的事實。
與上一個模擬相比，在模擬過程中原子的動作增大。如果你將虛擬熱浴溫度提升到更高的值（比如 3000 K），你將會發現晶體開始融化和部分非晶化。對晶體高溫退火實際上是一個產生非晶結構材料的通常方法。你將會在 創建無定形材料結構模型 這個例子中學到更多關於這個方法。
上面說過，調溫器（thermostat）時間標度參數可以控制向目標溫度接近過程的粒子速度。現在你將使用一個不同的調溫器時間標度值來重複同樣的模擬。打開上一個模擬的腳本，打開 MolecularDynamics 設置。保持其他設置不變，改變 Thermostat timescale 值為 500 fs。為了完全獲取這個持續很久的平衡相，增加步數到10000，改變文件名為 traj-nvt-nosehooverchain-1000K-500fs.hdf5。關閉設置，開始模擬。
模擬結果應該與下圖類似。
NVT Berendsen 和 Langevin 模擬
除了 Nose-Hoover chain 溫度器（thermostat）之外，你還可以選擇其他算法，比如 NVT Berendsen, Langevin, 或者NVT Nose-Hoover。最後的這個方法是 Nose-Hoover-Chain 方法更早和更簡單的版本，它的呈現更多的是出於歷史原因。如果可能，你應該優先選取 Nose-Hoover-Chain 方法，因為它給出許多更好的結果。Berendsen 調溫器（thermostat）提供一種算法有效的抑制溫度振蕩，使得在一些情況下系統更穩定。但是，相比於 Nose-Hoover-Chain，它不能精確得到正則系綜，儘管偏差非常小。由於這個原因，Berendsen 方法主要用在平衡態模擬，主要是在需要強健溫度控制的情況中。你可以為 NVT 模擬選擇的第三個調溫器（thermostat）類型是 Langevin 算法。這個算法求解 Langevin 方程，明確地在牛頓運動方程中包含摩擦以及隨機碰撞力以模擬粒子與虛擬熱浴之間的相互作用。不同於其他溫度器（thermostat）類型，Langevin 恆溫器（thermostat）將每個粒子分別與虛擬熱浴耦合。這產生一個很緊密的耦合，好像系統沉浸在一個虛擬的粘性介質中。但是它也以一種顯著的方式抑制了自然的動力學。出於此原因，Langevin 調溫器（thermostat）應主要被用來產生結構或者系綜取樣，而不是計算動力學性質。在 Langevin NVT 模擬中，耦合程度可以用 Friction 參數來設置。它被定義為量綱為時間的倒數，意味著它與 Berendsen 或者 Nose Hoover 溫度器（thermostat）的耦合參數的工作方式相反。增加 friction 值會導致與虛擬熱浴更緊密的耦合，和對動力學更嚴格的修正。
設置一個與 NVT 模擬包含相同模塊的腳本（BulkConfiguration, NewCalculator, Optimize ‣ MolecularDynamics）。
打開 MolecularDynamics 設置。通過查看 Type 下面的選項，你會找到兩個運行 NPT 模擬的選項：NPT Berendsen 和 NPT Melchionna 算法。
NPT Berendsen 方法的工作原理和 NVT 積分器相同，所以不能精確產生正確的物理系綜，儘管偏差通常也很小。而且，QuantumATK 中 Berendsen 的實現只考慮壓力標量 P=1/3(Pxx+Pyy+Pzz),（其中 Pij=−σij 為從應力張量 σ 得到的壓力張量的分量），視系統為各項同性。如果模擬像液體或者立方晶體這樣的各項同性系統，採用標量壓力是一個合適的選擇。對於各向異性系統，這個方法無法說明不同方向的不同模量。所以，NPT Berendsen 方法應只能被用於各項同性材料的平衡模擬。
為了將壓力分量 Pij與固定的壓力分別耦合，你更應該選擇使用 NPT Melchionna 方法7)。平衡定溫定壓 1 bar 下的矽超胞的合適的設置如下圖所示。
分子動力學（Molecular Dynamics）設置框包含了一般的設置，這些我們現在應該很熟悉。10000 個積分步長足以展示恆壓器（barostat）的效果。Log interval 可設置為 50 步以分析漲落。注意這將產生相當大的文件，或許你可以考慮在評估之後刪除軌跡文件。
另一個設置框包括恆溫器（thermostat）和恆壓器（barostat）的設置。恆溫器（thermostat）的設置我們已經在之前的部分討論過。設置External stress 為 1 bar。設置 Barostat time scale 為 100 fs。改變 Bulk modulus 為 98 GPa（矽的文獻值）。後面的參數不是特別重要，它只是告訴恆溫器（thermostat）體積如何隨著應力的變化而改變和確保壓力平衡的實際時間標度與所設時間標度的一致性。
最後，應力耦合編組框允許你選擇恆壓器（barostat）與哪個壓力/應力張量分量進行耦合，這意味著系統只與靜水壓力/應力耦合，而忽視剪切應力分量。這些設置對當前的例子十分有效。
要知道， QuantumATK NPT Melchionna 方法總是將所有壓力分量當做互相獨立的。並且，不僅是單個外界壓力，你也可以為每一個分量設定獨立的參考應力值來模擬力學剪切或者蠕變實驗。QuantumATK參考手冊 Reference Manual 和例子Simulating a creep experiment of polycrystalline copper 為此功能提供了更詳細的解釋。
將腳本送到 Editor 來增加一些額外的功能。為了可視化體積和壓力隨時間變化的漲落，將這幾行代碼附在腳本的末尾：
這幾行將從保存的軌跡（trajectory）中的每個快照中提取晶胞體積和應力，並畫出其值隨模擬時間的變化。
將腳本送到 Job Manager 來進行模擬。這將會花費幾分鐘來完成。使用 Movie Tool 中的 trajectory 打開軌跡，你就可以監測在模擬過程中晶胞的尺寸和形狀的變化。溫度在所選的模擬熱浴溫度300K周圍漲落。
你可以打開 Fluctuations_NPT_tau100.png 來查看體積的漲落。它應該與下圖類似。
上面的曲線顯示了晶胞體積顯著的振蕩。這個振蕩是 Nose-Hoover 類的溫度和壓力控制下的典型結果。只要他們的量級保持在合理的範圍內，並且在模擬過程中不會急劇增加，那麼就沒有必要擔心。振蕩的周期可以通過設置恆壓器（Barostat）時間尺度來改變。這個值選的越小，振蕩頻率越高。振蕩的振幅隨著恆壓器（barostat）時間尺度的減小而增加。
下面的曲線是模擬過程中的壓力，顯示了類似的行為。它在目標壓力 1 bar 附近振蕩，但是振幅較大。即時壓力值如此大的漲落其實在分子動力學模擬中很正常，因為壓力是基於應力來計算的。而應力對於體積變化非常敏感。此外，壓力是一個宏觀性質，對於如此小的系統考慮相應的平均值比即時值更合理。
NPT 模擬一個典型的應用是模擬相變，比如對石英在不同溫度下進行一系列模擬並測量平均晶胞體積隨溫度變化的函數，來研究α-石英到β-石英相變過程，如文獻8)所報導。
在進行分子動力學模擬過程中你可以固定選中原子的位置。設置這種約束最簡單的方式是打開 Molecular Dynamics 工具然後點擊 Add constraints 打開 Constraints editor 工具。
在工具欄裡你會找到一個表格列出了所有被標籤標記的原子組。你可以通過在展示窗口中（用滑鼠拖一個矩形或者點住 Ctrl 鍵選擇多個原子）選擇所需原子並點擊 Add tag from selection 來添加新的原子組。可以在 Constraint 欄中對每一個原子組指定約束條件。第一種將對應原子組裡的所有原子固定在它們的初始位置，而第二種將原子組內的原子整體視為一個剛性體。這意味著這些原子根據它們質心所受的力進行協同的移動。在當前的實行中只可能平移運動，而不能進行剛性體旋轉。
當使用約束時，你必須注意至少讓一個原子保持無約束的。像 NVT NoseHoover 和 NPT Melchionna 這樣的恆溫器（thermostat）需要至少兩個自由原子，因為它們還固定了系統的質心。恆溫器（thermostat）會自動適應由約束所導致的自由度數的減少。
然而，你需要知道 Movie Tool 或者 MD Analyzer 中一些功能不能與約束共同使用。比如，當分子動力學模擬中有約束時，Movie Tool 中顯示的溫度可能不會反應系統自由原子實際的溫度，而會顯示一個較低溫度。你可以通過在腳本裡添加動能（kinetic energy）代碼來監測溫度，如下例所示：
通常不建議在有約束原子的情況下實行定壓力/應力（比如 NPT）模擬，因為系統約束部分內產生的應力分布沒有從全局應力張量中移除。
在 DeviceConfiguration 中實現 NVE/NVT 分子動力學模擬本質上是和在 BulkConfiguration 中一樣的方式。拓展到中間區域的電極中的原子自動被約束不動以維持與電極的相容性。雖然技術上可以在 DeviceConfiguration 中實現定壓力/應力（NPT）模擬，由於約束的存在，不建議這種計算。對於 DeviceConfiguration 晶胞的優化有多種方法供選擇（參見例子 Geometry optimization of interfaces）。
1)
Tesoff: New empirical approach for the structure and energy of covalent systems. Phys. Rev. B 37, 6991, 1988
2)
W. C. Swope, H.C. Andersen, P.H. Berens, K.R. Wilson: A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 76, 637, 1982
3)
J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak: Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 81, 3684, 1982
4)
W. G. Hoover: Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695, 1985
5)
W. F. van Gunsteren, H. J. C. Berendsen: A leap-frog algorithm for stochastic dynamics. Mol. Sim. 1, 173, 1988
6)
G. J. Martyna, M. L. Klein, M. Tuckerman: Nosé–Hoover chains: The canonical ensemble via continuous dynamics . J. Chem. Phys. 97, 2635, 1992
7)
S. Melchionna, G. Ciccotti, B.L. Holian: Hoover NPT dynamics for systems varying in shape and size. Mol. Phys. 78, 533, 1993
8)
M. H. Müser, and K. Binder: Molecular dynamics study of the α-β transition in quartz: elastic properties, finite size effects, and hysteresis in the local structure. Phys. Chem. Minerals 28, 746, 2001