分子動力學并行算法的優化與應用

摘要：比較了小型機群中模擬小規模粒子系統采用原子分解和區域分解算法的模擬時間。針對原子分解算法，優化了各進程粒子的受力計算，尤其對嚴重影響并行效率的全局通信進行了優化，最終使得并行效率提高了30%以上。使用該優化算法對固態氬的導熱系數進行了計算，模擬結果與實驗值比較吻合，大大縮短了時間。

關鍵詞：分子動力學；并行算法；原子分解法；MPI_Allgather；鄰居交換算法；非阻塞通信

中圖分類號：TP3； TP301．6文獻標志碼：A

文章編號：1001－3695(2008)03－0718－03

分子動力學模擬[1]是目前在原子級對物體進行研究的主要方法之一。它將原子視為經典粒子，通過求解各粒子的運動方程得到不同時刻粒子的空間位置、運動狀態，從而統計出材料的宏觀行為特性。實驗研究表明，分子動力學模擬結果的精度主要取決于模擬系統所采用的分子間相互作用的勢能模型、模擬粒子規模以及時間步的多少。勢能模型主要取決于經驗；而粒子數和時間步越多，則精度越高，但同時增加了計算機的模擬時間。因此分子動力學并行算法的研究對非常耗機時的模擬具有重要意義。本文針對在實際應用中遇到的具體問題分析了適用于這些問題的分子動力學并行算法之一——原子分解法，并在此基礎上對其進行了優化。

1分子動力學并行算法

分子動力學模擬的步驟如下：a)設定模擬所采用的模型；b)給定初始條件；c)趨于平衡的計算過程；d)宏觀物理量的計算。在以上模擬過程中，c)即粒子的受力計算，是最費時的，一般占了整個模擬的90%左右，因此分子動力學的并行算法主要就是針對c)以實現計算時間的減少。

分子間相互作用力的計算花費了分子動力學模擬的大部分時間，其可分為長程力和短程力。前者是指所有粒子間相互作用，如電滲流的計算等；后者是指粒子間相互作用力隨距離增加而下降，在距離達到某個值后相互作用力可以忽略不計，該方法也是應用最廣的。假設分子間作用力為短程力，分子動力學并行算法主要分為三種[2]：a)原子分解法。將所有參加模擬的粒子編號平均分配給各計算節點處理器，它們只負責計算本進程擁有粒子。這種方法的好處是負載平衡容易保證、同步開銷小，但粒子在空間是隨機分布的，且在不停運動。為了使計算能夠進行，每運算完一個時間步后，需要將所有粒子位置、速度信息進行交換，這樣通信量會很大。b)分子間作用力分解法。將粒子受力計算的矩陣用子矩陣進行分割。這種算法既有原子分解法的優勢，又減少內存消耗和通信；但該算法編程復雜，不易作負載平衡，故一般不采用。c)區域分解法。將模擬系統劃分成多個物理子空間，每個計算節點對應一個(粒子分布均勻)或幾個(粒子分布不均勻有利于作負載平衡)子空間，這樣每個節點的粒子只需與周圍相鄰節點的粒子發生作用，只需要周圍鄰居的信息即可，所以通信量比較小。但由于每個區域的粒子數不同，負載不容易平衡，同步開銷較大。該方法在大規模系統中應用最為廣泛。

2原子分解算法的優化

本文實現的分子動力學并行算法中采用的是原子分解法而非應用較為廣泛的區域分解法，主要原因基于以下幾點：a)串行程序并行化使用原子分解法最為簡單方便。b)許多分子動力學模擬系統一般在數千以上規模，但是由于時間步長在飛秒級，模擬步數一般都是幾百萬，所耗機時長時需以周來計算。區域分解法主要適用于大規模粒子系統，而且是在大型機群系統中。對于在小型機群中的小規模系統的模擬采用原子分解法可以得到比較理想的并行效率[3]。區域分解法并不具備優勢，從第3章給出的小型機群中區域分解與原子分解法時間比較可以看出。c）一些模擬必須考慮分子間的長程作用力，如電滲流的計算，而其他兩種并行算法均不適合，只能采用原子分解法。d)許多模擬中勢函數使用了三體、四體作用勢，使用區域分解不但復雜，效率也會很低，在此類問題的處理上許多大型分子動力學模擬軟件，如CHARMM和GROMOS也是采用原子分解的方法[2]。因此編寫出基于原子分解法性能優良的并行分子動力學代碼對于存在以上幾種情況的分子動力學研究具有重要意義。

對于模擬同樣的系統，幾種算法的計算時間都是一樣的，但在通信時間上原子分解、力分解、區域分解依次減少，所以原子分解法的效率主要取決于全局通信在通信與計算總時間中的比例。通信時間所占比例越小，則其并行效率越高。因此減少通信時間和增加計算量是提高效率的兩種途徑。在計算量一定且相對較小時，通信對并行效率的影響較大，這時應主要考慮在全局通信上進行優化。此外，在計算上也可以作一些優化。在本文的串行程序中，對于粒子間作用力的計算采用牛頓第三定律，這樣兩粒子間作用力只要計算一次就可以了，節約了一半的時間。但是在采用并行算法時，如果采用牛頓第三定律，通信則會變成原來的兩倍，對機群系統是不可取的，而且在力的模型不是很復雜時效果也不是很理想[2]，所以在并行算法中沒有采用牛頓第三定律。這樣，各進程N/P(N是粒子總數，P是進程數)個粒子雖然增加了cache命中率，但是最終計算效率的提高也只是略大于50%。如果將各進程粒子與本進程和其他進程送過來的粒子計算分開，前者采用牛頓第三定律，這樣計算效率可以提高60%左右，效果比較理想。而且對于模擬對象固體的情況，由于粒子運動基本是在原位置作小幅振蕩，可以采用將粒子分解與區域分解結合的方法，在分配粒子時按照其空間位置進行分配[4]，以增加本進程內的粒子計算量，減少與其他進程內粒子的計算，從而提高計算效率。對于液體，由于粒子運動較為劇烈，可以間隔n步對粒子進行一次重新分配，n可以根據測得的計算時間的顯著變化來確定或對各進程粒子進行不定時監測，根據運動出本進程區域的粒子數量來確定n。然而盡管如此，再加上全局通信時間，最終的并行效率也只是在60%左右。此時只能考慮在全局通信上作進一步優化。

本文采用的并行編程環境是消息傳遞并行編程環境MPI[5，6]。一般情況下，當各進程粒子數相等時，在程序中調用MPI中的組通信函數MPI_Allgather來實現所有粒子信息的收集，不等則調用MPI_Allgatherv。本文僅討論各進程粒子數相等的情況，為下一步的計算提供數據。在程序中直接調用MPI_Allgather可以方便編程，提高程序的可讀性、移植性以及通信效率。該語句也是MPI庫中使用頻率最高的集合通信函數之一。最新的MPICH版本使用了三種算法實現MPI_Allgather，即環、遞歸倍增和Bruck[7]。在調用MPI_Allgather函數時，MPI系統根據消息的長度、參與進程的數量、網絡帶寬及網絡延遲選擇最優的算法。此外，針對Allgather算法，國內外許多學者對此進行了研究和優化[7，8]。陳靖等人[9]提出了一種鄰居交換算法，在本文中也被用來與其他三種算法一起在具體機群上進行測試以確定在本文中使用何種最優算法。

鄰居交換算法因其是相鄰進程兩兩交換消息，具有良好的通信局部性，可以最小化TCP/IP通信量，能夠獲得較好的通信性能。對于P（偶數）個進程參與時，該算法需要P/2步通信完成全收集的工作：a）偶數號進程給其右邊鄰居((i+1)% p)發送自己擁有的一個數據塊；奇數號進程則向其左邊鄰居((p+i+1)% p)發送一個數據塊。b）反之，偶數號進程與其左鄰居交換消息，共有自身的和a)接收到的兩個數據塊。c）后邊奇數步重復a)的操作，只是交換的數據都是上一步接收到的兩個數據塊，偶數步則同b)。

對全局通信作優化，也就是對MPI_Allgather進行優化，是在針對具體模擬系統和測試環境得到的最優算法基礎上，使用點對點通信函數編寫的代碼實現該算法以取代直接調用MPI_Allgather組通信函數，對于該點對點通信使用非阻塞通信取代MPI_Allgather調用的阻塞通信算法。MPI_Allgather幾種算法的通信采用的都是阻塞通信。這意味著在每一步通信完成前，CPU只能處于等待狀態，消耗了大量機時，尤其對于通信時間在總時間中占較大比重的情況。在實際測試過程中，比例平均在30%以上。如果改用非阻塞通信，則CPU在每一次通信的開始提出通信請求，然后立即返回，通信交給特殊硬件去完成，這時可以繼續執行其他指令。可以在通信的同時進行與通信數據無關的計算任務，實現計算與通信的重疊，大大縮短總的時間。圖1是采用原子分解法，阻塞與非阻塞通信對比的流程圖。

各進程粒子與本進程內粒子的作用力的計算與需通信獲得的數據無關，可以在非阻塞通信時進行計算。如果通信時間長于本進程內粒子的計算，還可以增加與通信所獲數據無關的其他計算，即各進程粒子與上次通信獲得的粒子的作用力的計算可以與后邊通信實現重疊。本文中僅使用了通信與本地進程粒子計算的重疊。每步計算的量根據每步通信時間來分配，使得計算時間與通信時間相差盡可能小，以實現最大程度的重疊。

3測試結果

為驗證該優化算法的正確性及有效性，筆者在東南大學計算機學院高性能計算中心對其進行了測試。該機群包括32個計算節點，其主要計算節點的配置為雙Xeon 3.0 GHz CPU，4 GB內存，互聯網絡為10 GBps Infibiband。

在該機群上，對于中等長度8 KB左右的消息、全收集MPI_Allgather的四種算法中，以鄰居交換算法的性能最優，最終也是采用這種算法對程序進行了優化。首先，比較了八個進程時使用優化后原子分解算法和區域分解算法對于小規模粒子系統的模擬時間。本文編寫的區域分解算法代碼只考慮了八個進程以上的情況，并且對通信作了與計算重疊的優化，對固態氬的受力進行了模擬。比較結果如表1所示。可以看出，模擬10 000以內小規模粒子系統時使用原子分解算法是優于區域分解算法的。

然后對原子分解算法進行了具體測試。分別用1、2、4和8個進程對擁有6 912個粒子的固態氬的分子動力學模擬時間進行了測試，得到了串行、優化前并行和優化后并行程序在每個時間步所耗費的時間。從圖2的單步時間比較可以看出，優化后并行效率得到了大幅提升，同時隨著進程數的增加，并行效率有所降低。這主要是因為進程數增加，本地進程的粒子數減少，計算效率下降，而且通信時間有所增加。該算法在小型機群系統中可以獲得比較理想的并行效率。

采用四個進程分別對不同小規模的粒子系統模擬所耗費的時間進行測試，與串行程序的模擬結果進行比較得到了并行效率，如圖3所示。可以看出，隨著粒子數的增加，并行效率不斷增加，這是因為粒子的增加提高了計算效率。盡管粒子數增加的同時通信時間也在增加，但是增加的通信在使用非阻塞通信時可以獲得與計算更好的重疊。

為了驗證其正確性，用該算法對固態氬的導熱系數采用平衡態分子動力學進行了計算，氬原子間相互作用力采用LennardJones作用勢[1]，勢能參數為σAr=3.4 A。，εAr/kB=121．0 K。對于4 000個粒子的體系，在17 K的溫度下模擬50萬步，前20萬步是系統的平衡時間，最終得到的結果與實驗值[10]吻合較好，如圖4所示；模擬時間也大大縮短，串行程序時間是3.9 h。而對于優化后并行程序在使用四個進程時總的模擬時間是1.25 h，并行效率可達78%。

4結束語

本文通過測試比較得出在小型機群中模擬小規模粒子系統原子分解算法要優于區域分解算法。在此基礎上，對原子分解法進行了分析和優化，分析了其應用場合及造成效率低的主要原因；針對這些原因進行了計算和通信方面的優化，主要針對嚴重影響并行效率的全局通信，使用非阻塞通信實現通信與計算的重疊。根據全局通信幾種算法的介紹[9]，使用筆者編寫的點對點通信代碼替代直接調用MPI_Allgather語句，在所測試的機群系統中測得使用鄰居交換算法性能最好，最終選擇了該算法。在進行非阻塞通信的同時計算各進程粒子與本地進程粒子間及上次通信獲得的其他進程粒子間的相互作用，獲得了比較好的效果，并且通過計算固態氬的導熱系數證明了該優化程序的正確性和有效性。

致謝:本文的測試工作得到了東南大學計算機科學與工程學院高性能計算中心的支持，在此表示衷心感謝。

參考文獻：

［1］ALLEN M P，TILDESLEY D J.Computer simulation of liquids[M].Oxford: Clarendon Press，1987.

[2]PLIMPTON S.Fast parallel algorithms for shortrange molecular dynamics[J].Journal of Computational Physics，1995，117(1): 119.

[3]吳江濤，劉志剛，董小社，等.RS/6000機群系統中分子動力學并行算法的研究[J].工程熱物理學報，2002，23(4): 418－420.

[4]BROWN D，CLARKE J H R，OKUDA M，et al.A domain decomposition parallelization strategy for molecular dynamics simulations on distributed memory machines[J].Comput Phy Commun，1993，74(1):67－80.

[5]都志輝.高性能計算之并行編程技術——MPI并行程序設計[M].北京：清華大學出版社，2001.

[6]莫則堯，袁國興.消息傳遞并行編程環境MPI[M] . 北京：科學出版社，2001.

[7]BRUCK J，HO C，KIPNIS S，et al.Efficient algorithms for alltoall communications in multiport message passing systems[J].IEEE Trans on Parallel and Distributed Systems，1997，8(11): 11431156.

[8]RAJEEV T，WILLIAM G.Improving the performance of MPI collective communication on switched networks[C]//Proc of the 10th European PVM/MPI User’s Group Meeting.Venice:[s.n.]，2003:257－267.

[9]陳靖，張云泉，張林波，等.一種新的MPI_Allgather算法及其在萬億次機群系統上的實現與性能分析[J].計算機學報，2006，29(5): 808－814.

[10]CHRISTENT D K，POLLACK G L.Thermal conductivity of solid argon[J].Phys Rev B，1975，12（8）:3380－3391.

“本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文”