分子動力學模擬研究的并行計算程序的改造

分子動力學模擬研究的并行計算程序的改造

1物理模擬測試在材料微觀結(jié)構(gòu)的研究中，由于實驗條件的限制，很難獲得所需的結(jié)果。例如,對液態(tài)或非晶態(tài)結(jié)構(gòu)的短程序,就難以用X射線衍射反映出其確切的微觀結(jié)構(gòu);在由液態(tài)金屬急冷到非晶態(tài)的凝固過程中,就根本無法用儀器觀測到其微觀結(jié)構(gòu)的瞬時變化等等。然而,隨著計算機技術(shù)的飛速發(fā)展,用計算機模擬研究的方法,卻往往能夠比較方便地得到這些結(jié)果,因而它已日益成為材料微觀結(jié)構(gòu)研究的一種重要手段。計算機模擬有多種方法,其中分子動力學的方法比較適合于液態(tài)及非晶態(tài)金屬動態(tài)過程的研究(我們早在九十年代初就開始進行這方面的研究工作),至今仍是液態(tài)及非晶態(tài)金屬微觀結(jié)構(gòu)研究的一種重要方法。當這種方法在中型機或微機上進行時,由于其內(nèi)存和運算速度的限制,模擬研究只能限于500～1000個金屬原子的小系統(tǒng)。因而模擬得到的結(jié)果雖然也能揭示一些微觀結(jié)構(gòu)的特征和規(guī)律,但與實際的大系統(tǒng)情況有較大的差異。如果我們能夠在并行處理系統(tǒng)上對更大量的原子系統(tǒng)進行模擬研究,其結(jié)果必然會更接近于實際,從而對生產(chǎn)實踐將會更有實際指導意義。然而在并行處理系統(tǒng)上進行大型模擬計算必須采用并行計算才能實現(xiàn),才能充分發(fā)揮并行處理系統(tǒng)的優(yōu)勢。因此,本文在原有研究工作的基礎(chǔ)上,嘗試對原有分子動力學模擬計算程序進行并行化處理,使所研究的系統(tǒng)的原子數(shù)達到50000～100000個以上,從而使模擬研究所得的結(jié)果與真實情況更為接近。2原子數(shù)計算的一般程序分子動力學是求解相互作用粒子體系的經(jīng)典運動方程的理論。只要知道所研究系統(tǒng)中原子的初始位置及速度,以及原子間的相互作用勢V,就可以計算出每個原子所受的力,從而計算出每個原子以后的運動,得到原子的新位置及速度。只要一步步地跟蹤計算下去,就可以模擬出系統(tǒng)中原子運動的情況,得到系統(tǒng)瞬時的微觀結(jié)構(gòu)組態(tài),并可由此計算出該原子系統(tǒng)的熱力學統(tǒng)計性質(zhì)。在分子動力學模擬中,用計算機對粒子體系的經(jīng)典運動方程求解時,主要是由VERLET算法實現(xiàn)的。由于VERLET算法比較成熟,應(yīng)用較廣,并且我們已由此得到了許多少量原子(500～1000個)系統(tǒng)的模擬結(jié)果,因此本文在考慮并行算法時,主要是在原算法的基礎(chǔ)上加入并行部分,而不是完全從并行的角度重新設(shè)計算法,要求所得到的結(jié)果(50000個原子系統(tǒng)的)能與以前的結(jié)果有可比性。在VERLET算法中,如果已知t時刻及前一步t-dt時刻原子的坐標r(t)和加速度a(t),則在t+dt時刻,將有r(t+dt)=2r(t)-r(t-dt)+dt2a(t)(1)速度可以不計算,但若要估算動能,則速度可由下式得到:v(t)=[r(t+dt)-r(t-dt)]/(2dt)(2)其中加速度a的計算即是原子受力f的計算。計算f采用雙體作用勢近似,即第i個原子受到的力為fi=Σfji(t)?i≠j(3)fi=Σfji(t)?i≠j(3)其中fji(t)是第j個原子對第i個原子的作用力。由(3)式容易看出,從理論上來說,其中f的計算量是與原子數(shù)的平方成正比,而其它量的計算量一般是與原子數(shù)成正比,因而f的計算將占用計算中的大部分時間,并且比例將隨原子數(shù)的增加而增加。若系統(tǒng)取500個原子,一般計算力的程序占時為85%以上?？梢钥隙?隨著系統(tǒng)選取的原子數(shù)的增加,其比例將會更大。在原程序的實際計算中,采用了一種犧牲空間換取時間的作法,即在t時刻記錄下每個原子周圍RCUT距離內(nèi)的原子(超過此距離的原子則認為對其影響可忽略,即其雙體勢將趨近于0),由近鄰原子表數(shù)組記錄,而這些原子簡稱為近鄰原子。如果在下一步t+dt時,各原子的位移均未超過某一限定值,則認為各原子周圍的近鄰原子仍是t時刻記錄的那些近鄰原子,此時程序調(diào)用FORCE2計算各原子的受力;若上述情況不滿足,則某一原子受力的計算需要考慮其它所有原子,此時程序?qū)⒄{(diào)用子程序FORCE1來計算力,同時各原子周圍的近鄰原子發(fā)生變化,而要求更新近鄰原子表數(shù)組。由于與某個原子的距離在RCUT內(nèi)的原子數(shù)一般不會隨整個系統(tǒng)原子數(shù)的增加而發(fā)生較大的變化,因此FORCE2的實際計算量將與系統(tǒng)的原子總數(shù)成正比。而FORCE1要求計算兩兩原子之間的作用力,其計算量將會與系統(tǒng)的原子總數(shù)的平方成正比。實際上整個計算力的計算量M=λ1M1+λ2M2,其中M1、M2分別是FORCE1和FORCE2的計算量,而λ1、λ2則分別是FORCE1和FORCE2被調(diào)用的次數(shù)。當原子總數(shù)增加或者dt增加時,FORCE1的調(diào)用次數(shù)將越來越多,而FORCE2則相反,甚至當dt超過某一值時,FORCE2將不會被調(diào)用。采用這種策略將節(jié)省計算時間,但須記錄下每個原子周圍RCUT內(nèi)的原子情況(比較多時會達500～600個),將會占用很大的內(nèi)存,在原子總數(shù)很大的情況下,內(nèi)存將難以滿足需求,這也是制約大量原子系統(tǒng)模擬計算的一個重要因素。基于以上分析,本文采用以下的并行算法:將計算力的部分程序作并行化計算處理。采用MASTER-SLAVE方式,將所有原子的數(shù)據(jù)按原子編號平均分到各SLAVE結(jié)點機,包括每個原子的坐標、近鄰原子表數(shù)組(即某原子周圍RCUT范圍內(nèi)的原子列表)等數(shù)據(jù)。例如,整個系統(tǒng)有10000個原子,分到10臺處理機上,則1～1000號原子分到1號機上,而1001～2000號原子分到2號機上,如此類推。當主程序要計算力時,將所有原子的當前坐標廣播到各個處理機,然后由各個處理機同時計算其所屬范圍內(nèi)的原子的受力。各個SLAVE計算完以后,再將結(jié)果(即各原子的受力)發(fā)回到MASTER,由主程序繼續(xù)下一步。同時,各原子的近鄰原子表分開存放在各SLAVE機中,即使有所變化,也只由各SLAVE機完成,而與MASTER無關(guān),這樣主程序中占用龐大內(nèi)存的近鄰原子表數(shù)組可以去掉。采用上述算法有以下幾個優(yōu)點:(1)算法未對原算法進行本質(zhì)性的修改,只是加入了并行部分,所得結(jié)果與原來一致,便于檢驗和比較計算結(jié)果;(2)算法采用MASTER-SLAVE方式,簡潔明了,易于維護,可隨時調(diào)整并行機數(shù)量;(3)算法解決了原算法中的內(nèi)存問題,將原來單機上需要的大量內(nèi)存分開到各個并行機中,從而使可計算原子數(shù)提高了很多。上述算法的理論加速比是R1*C/2,C是SLAVE機數(shù),R1是一個系數(shù),指計算力的部分在整個計算中所占的比例,當原子總數(shù)很大時,R1一般在90%以上。另外再除以2是因為單機計算時,兩原子之間的相互作用力是相對的,只須計算一次,而并行計算時,由于分開計算,將會重復(fù)計算一次。3sla智能中的并行化處理本程序在某并行處理系統(tǒng)上運行,采用MASTER-SLAVE方式的PVM編程實現(xiàn),將原單機程序中計算力的部分,即子程序FORCE1,FORCE2中的依次計算各個原子受力的部分(為主要部分)作并行化處理,轉(zhuǎn)化成SLAVE程序。每次計算力時,主程序調(diào)用子程序FORCE1或FORCE2,而子程序再依情況將所需要的數(shù)據(jù)傳給各個SLAVE。各個SLAVE收到數(shù)據(jù)后即開始這一輪計算,結(jié)束后將結(jié)果回傳給MASTER,然后處于等待狀態(tài)中,直到下一輪計算力的部分開始。MASTER依次接收各個SLAVE的數(shù)據(jù),然后繼續(xù)下一輪,直到整個計算結(jié)束。4al原子的密度分析結(jié)果經(jīng)過并行處理后的計算程序已在某并行處理系統(tǒng)上穩(wěn)定運行。我們使用50個處理單元進行并行計算,實現(xiàn)了50000～100000個金屬原子的大系統(tǒng)的模擬計算。這在國內(nèi)材料科學領(lǐng)域的分子動力學模擬中尚屬首次,在國際上亦處于先進水平?，F(xiàn)已計算了50000個Al原子在943K時等溫平衡的微觀結(jié)構(gòu)情況,簡要介紹如下:(1)從圖1所示的雙體相關(guān)分布函數(shù)g(r)的曲線可以看出:模擬計算所得曲線與實驗曲線符合較好。由于g(r)是一個統(tǒng)計平均的結(jié)果,所以500個原子的g(r)與50000個原子的相差甚微。而正是這種統(tǒng)計平均,使得許多微觀結(jié)構(gòu)的重要差別被淹沒掉,為此必須作進一步的深入研究。(2)再從兩者的密度分析結(jié)果來看(如表1所示),50000個原子的漲落比500個原子的要大,這是一個更加符合實際的結(jié)果。由此可向我們提供更多的微觀結(jié)構(gòu)信息。如50000個原子系統(tǒng)中可以找到含幾百個原子的團簇結(jié)構(gòu),這在500個原子的小系統(tǒng)中顯然無法找到。5并行算法的實現(xiàn)隨著計算機技術(shù)及分子動力學模擬技術(shù)的快速發(fā)展,在物理、化學、生物及材料科學等領(lǐng)域,對大量原子、分子系統(tǒng)進行計算機模擬研究的熱潮必將會出現(xiàn)。而從目前來看,模擬計算只有采用并行算法才可以實現(xiàn)。本文所述的算法已經(jīng)能夠得到50000個原子系統(tǒng)的模擬結(jié)果,但只是做出了這方面的初步嘗試