分子動力學模擬算法框架

技術背景

分子動力學模擬在新材料和醫葯行業有非常重要的應用，這得益於分子動力學模擬本身的直觀表述，用宏觀的牛頓力學，結合部分微觀的量子力學效應，就能夠得到很好的符合統計力學推斷的結果。可以說，分子動力學模擬從理論上跨越了物理化學生物等多個學科，而從實踐上又包含了計算機科學、人工智能的大量輔助和優化，綜合性非常強。越是綜合性強的研究方向，就越有必要梳理清楚其主干脈絡和工作流程。

1 坐標和速度的初始化

當我們需要研究一個特定的體系時，首先需要指定這個體系的初始空間和坐標，再進行演化。初始坐標的獲取方式有很多種，比如從數據庫中直接獲取、根據拓撲結構用特定的工具進行生成（比如之前的博客介紹過的OpenBabel），還有就是在材料學領域常用的根據晶體結構去生成有規律的體系初始坐標。初始速度的選取方式，其實可以選擇全0的初始狀態，讓體系在力場作用下進行自由演化；也可以給全隨機的初始速度，讓一開始的體系能量和溫度達到一個比較高的值，有可能可以提高模擬的效率；還有一種在材料學領域比較常用的方法，就是指定一個初始溫度，通過能均分定理對體系中的各個原子進行速度分配。如下所示是一個通過OpenBabel生成一個水分子初始坐標的簡單示例：

$ obabel -:O --gen3d -omol2 -O H2O.mol2
1 molecule converted
$ cat H2O.mol2
@<TRIPOS>MOLECULE
*****
 3 2 0 0 0
SMALL
GASTEIGER

@<TRIPOS>ATOM
      1 O           0.9606    0.0592    0.0649 O.3     1  HOH1       -0.4105
      2 H           1.9285    0.0383    0.1064 H       0  HOH0        0.2052
      3 H           0.6818   -0.3573    0.8942 H       0  HOH0        0.2052
@<TRIPOS>BOND
     1     1     2    1
     2     1     3    1

應該說，在坐標的初始化中，由於分子本身的性質很大程度上是由其拓撲連接性所決定的，所以生成分子初始坐標要以拓撲結構為依據。而初始速度的選取，更多的是一種經驗方法，只要初始速度值不要過大（否則會使體系劇烈震盪甚至奔潰），怎么選取初始速度就只是影響最終的模擬效率。比較典型的一個案例就是通過Shooting的方法對化學反應過程進行加速，其實也是對初始速度的一個調整。

2 Neighbor List與Cell List

在計算非鍵相互作用時，如果對一個\(N\)原子的體系進行遍歷，可達\(\frac{N(N-1)}{2}\)的計算量。由於模擬的體系一般都較大，\(O(N^2)\)的復雜度屬實有點高了。因此為了簡化計算量，一方面通過球形截斷來選取近鄰表（如下圖所示，圖片來自於參考鏈接1）：

僅在截斷的球形范圍內考慮兩兩相互作用。再對空間進行切割，構建Cell List，只要設定好Cell List的大小，比如設定正方體Cell List的邊長等於截斷半徑，這樣一來，我們在統計一個Cell中的近鄰關系時，只需要關注總共周邊的27個Cell中的原子即可，大大縮減了計算量。Cell List示意圖如下所示（圖片來自於參考鏈接2）：

關於Cell List的實現，此前筆者專門寫過一篇博客用於介紹在Python中的GPU加速打格點算法實現，大致效果如下：

由於是用了GPU來計算，因此有較好的性能表現。一般情況下，構建近鄰表的這一步是對GPU的顯存要求最高的步驟之一。即使是使用了Cell List的方法，一般也不會每一步都去更新近鄰表。常用的策略是，設置一個interval，比如5-10，也就是每5-10步更新一次近鄰表。而且也不需要再去檢索整個鄰近的Cell中的所有原子，只需要監測skin region中的原子變動即可。這也是設立一個緩沖區skin region的原因，當然，我們需要保障在指定的interval下原子不會跑的太遠。

3 獲取作用力

獲取作用在每一個原子或者團簇上的作用力的方法有很多種，比如量子化學的ab-initio的方法、擬合力場的方法以及近幾年發展出來的深度學習的方法，都有非常多的應用且各有優勢。從頭算的量子化學方法，精度非常高，但是耗時較長。擬合力場速度非常快，精度也在可接受的范圍內，是傳統分子動力學模擬中應用最廣泛的方法。但是這兩個方法都要面臨一個問題：需要手動微分，要么就是復雜的解析推導，要么就是數值差分。而現在計算機領域的各種深度學習框架，都已經提出了最新的解決方案：自動微分。因此，使用深度學習的方案去構建力場，在充分的訓練集和合理的損失函數下，可以達到非常快的速度以及可以接受的精度。還有兩個非常重要的點：整個過程只要使用深度學習的算子，就是完全連續可微的，而且往往這些深度學習框架都對硬件系統有非常好的可拓展性和可遷移性。

除了上述的三種方法之外，結合這幾年非常火的量子計算，又出現了新的基於量子計算機的分子動力學模擬框架，比如參考鏈接3和參考鏈接4。目前來說已經提出來的主要框架還是基於Variational Quantum Eigensolver變分量子求解器（VQE），通過VQE可以計算出來類似於量子化學從頭算方法的勢能面，進而根據費曼-赫爾曼定理做差分得到原子上的作用力，就可以用於分子動力學的模擬了。由於這里細節較多，在本文不進行展開講解，需要后面專門開一篇文章介紹VQE的框架。簡單點評一下這兩個基於VQE的分子動力學模擬的工作的話，就是只適用於當下量子芯片的Demo使用，距離真在在分子動力學模擬場景下的應用還非常的遙遠。而VQE框架本身是NISQ近期量子計算時期下所誕生的產物，其最promising的點在於，可以在近期含噪的量子芯片下實現一個比較大規模體系的分子的勢能面的計算（計算量已經超越了傳統計算機的范疇），但是實際操作中在精度和速度上是否能有優勢，暫時還是存疑的狀態。另外，基於類似的變分量子計算的框架，還出現了量子機器學習這個兩大buff加持的算法框架，同樣可以參考傳統深度學習的方法實現一個量子計算版本的神經網絡來替代力場，目前似乎還沒人這么做。

4 積分器/體系演化

在具備了初始狀態的坐標和速度，以及上一步計算得到的原子作用力之后，就可以使用牛頓力學讓體系向前演化一步。其中涉及到的算法內容，可以參考上一篇博客從劉維爾方程到Velocity-Verlet算法，其中主要就是Verlet、LeapFrog和Velocity-Verlet三個算法，在分子動力學模擬中也叫積分器。這個部分的計算量並不是很大，但是會導致體系的運動，如果需要保存中間過程，注意每隔一定的模擬步數，要把軌跡保存下來。

5 判斷控制/約束條件

上一個章節中說的是體系的運動，那么在一定的時間步長（通常是一個飛秒）下，原子運動的幅度大小，是跟多個因素相關的：1. 原子本身收到的作用力大小；2. 原子的質量；3. 在控溫控壓，有可能受到系統粘滯作用力的影響或者碰撞影響，也有可能被縮放；4. 對於加了約束條件的過程，還會受到施加的約束作用，比如等效的鍵向作用力（施加約束的過程，本身就是為了控制部分原子之間相對位置大小）。因此，我們在分子動力學模擬的演化過程中，還需要考慮這些控制和約束的方法。

常用的控溫技術有：Andersen Thermostat速度重置方法、Berendson Thermostat速度縮放方法和Extended Ensemble擴展系綜方法。但是Andersen Thermostat存在相空間軌跡不連續的問題，Berendson Thermostat可能會導致不合理的相的出現，因此最常用的就是幾種擴展系綜的方法：Nose-Hoover Thermostat、Nose-Hoover Chain Thermostat和Langevin Dynamics Thermostat，這些方法的基本原理就是在系統哈密頓量上加上溫度的相關項，這樣就可以驅動系統演化到一個溫度恆定的本征態上，以達到控溫的目的。常用的控壓技術有：Volume Rescaling體積重設方法、Berendsen Barostat體積坐標縮放方法和Extended Ensemble。類似於控溫技術的，最常用的還是在哈密頓量上加一個體積項，通過驅動系統演化到一個穩定體積所對應的本征態，從而達到控制體積的目的。而這里提到的約束條件，相對的就更加具體，常見的比如LINCS算法、SETTLE算法、RATTLE算法和SHAKE等，其基本思想就是固定某兩個原子之間的距離保持不變。用一種比較形象的說法就是，在兩個原子之間加了一根沒有質量的柱子限制兩者的相對運動。如下所示是SETTLE算法的實現，先不加約束的位移一步，然后再把約束加進去，同時對坐標和速度進行更新：

6 結果輸出與增強采樣

經過上面的步驟，我們其實已經完成了一步的分子動力學模擬，如果此時達到了結果輸出的點，應該將坐標軌跡保存到一個文件中。需要特別說明的是，有一些格式並不支持持續的數據寫入，比如npz格式，每寫入一次，就會自動把前面寫入的內容全都覆蓋，因此要挑選好存儲格式，常用的是二進制格式的文件或者文本格式的文件。除了存儲完整的坐標軌跡之外，還需要存儲一些必要的參數，比如單步的速度等，這樣即使程序終止，下一次啟動程序的時候也可以從斷點開始，而不需要重新開始。

而增強采樣（Enhanced Sampling）也是一個基於結果來分析的工作，稱之為后處理。雖然是后處理，但是這個處理的結果會反饋到下一步的迭代中去，作為一個獨立的等效勢能去推動原子的運動。常見的工具是PLUMED，前面寫過一篇相關的安裝與使用的文章。當然其實采樣算法也並不是難以實現，更大的難點是在於跟已有框架的結合與交互，相對來說門檻要高一些。常用的增強采樣的方法比如MetaDynamics，可以用這樣一個比較具象的案例來描述：在一個自定義的序參量所對應的勢能面上堆土，直到分子體系可以在給定的本征態之間快速切換（勢能差距較低）。大致結果如下圖所示（圖片來自於參考鏈接5）：

在通過新的坐標更新了力場參數之后，加上增強采樣所對應的等效勢能產生的等效作用力，就可以開始下一輪的迭代，就重新回到了第3步，不斷的迭代，直到達到目標的迭代步數。一般情況下分子動力學沒有特別明顯的停止迭代的指標，更多的是設定一個時間閾值，達到閾值就停止迭代更新。

總結概要

分子動力學模擬是一個跨越眾多學科領域的強大工具，從物理學的角度來看分子動力學模擬的話，其基於量子力學（量子化學）構建模型，通過牛頓力學進行演化迭代，最后能夠在時間平均上等同於統計力學的系綜平均，是一個堪比復變函數歐拉公式的優美過程。本文就當前分子動力學模擬的框架進行了整體介紹，其中並不展開講解各項技術內容，但是也為感興趣的研究人員提供一個簡單的入口。

版權聲明

本文首發鏈接為：https://www.cnblogs.com/dechinphy/p/md-framework.html

作者ID：DechinPhy

更多原著文章請參考：https://www.cnblogs.com/dechinphy/

打賞專用鏈接：https://www.cnblogs.com/dechinphy/gallery/image/379634.html

騰訊雲專欄同步：https://cloud.tencent.com/developer/column/91958

參考鏈接

1. http://blog.sina.com.cn/s/blog_78f0b32101011mlj.html#:~:text=近鄰列表：近鄰列表 (neighbor,list，也叫做Verlet list)是分子動力學模擬中為了加速體系粒子間相互作用的力的計算而設置的一種數據結構。
2. https://www.jsjkx.com/CN/article/openArticlePDF.jsp?id=512
3. https://aip.scitation.org/doi/10.1063/5.0046930
4. https://pubs.acs.org/action/showCitFormats?doi=10.1021/acs.accounts.1c00514&ref=pdf
5. https://zhuanlan.zhihu.com/p/95824483