基于變步長Velocity Verlet 算法的盧瑟福散射模擬的優(yōu)化

周 偉，張宜新，劉學(xué)鋒，楊喜峰，王殿生，王玉斗

（中國石油大學(xué)（華東） 理學(xué)院，山東 青島 266580）

盧瑟福散射是近代物理科學(xué)發(fā)展史上一個非常重要的實驗。該實驗驗證了原子的核式模型，奠定了原子物理學(xué)的基礎(chǔ)[1-2]。同時，盧瑟福散射還推動了許多技術(shù)應(yīng)用方面的進步。例如：利用離子束與物質(zhì)相互作用對物質(zhì)進行元素成分分析和結(jié)構(gòu)分析的離子束分析技術(shù)[3]；利用盧瑟福散射中的大角度偏轉(zhuǎn)進行物質(zhì)分析的背散射技術(shù)[4]等。

由于盧瑟福散射有著很高的技術(shù)價值，因此它一直保持著較高的研究熱度。計算機模擬是一種常用的研究盧瑟福散射的方法。目前，國內(nèi)外有許多關(guān)于計算機模擬盧瑟福散射的研究。文[5—6]中，作者直接應(yīng)用盧瑟福散射公式計算了單靶核對α 粒子的散射幾率。然而這種方法無法模擬出α 粒子的運動情況，更無法研究實際散射過程中存在的多靶核的影響。文[7]中，作者利用Visual Basic 語言完成了盧瑟福散射軌跡的模擬，模擬出了不同α 粒子入射能量和瞄準距離下的α 粒子散射的運動軌跡。但其主要內(nèi)容是單個粒子散射軌跡的演示，并沒有討論散射幾率。文[8—9]中，作者利用歐拉方法與蒙特卡羅方法模擬了散射過程，既模擬了散射的運動軌跡，又得到了散射概率曲線。但是，當(dāng)靶核數(shù)目較多或模擬次數(shù)較大時，文中所使用的程序出現(xiàn)了效率低、計算時間長的問題，難以滿足研究要求。

目前，盧瑟福散射的計算機模擬效率低主要有兩方面原因：一是核心算法的精度偏低，導(dǎo)致了單位步長較??；二是模擬中沒有考慮不同條件下粒子受力的變化，采用了單位步長固定的模擬方法。鑒于以上不足，本研究基于文[9]中所使用的模擬框架，采用了新的優(yōu)化算法，并在此基礎(chǔ)上引入了動態(tài)步長，大幅度提高了盧瑟福散射的模擬效率。

1 模擬方法

首先介紹本文中進行盧瑟福散射的計算機模擬的基本流程，然后詳細介紹幾種軌跡模擬算法的原理并比較其效率，最后給出提高模擬效率的變步長處理方案。

1.1 模擬過程

本研究將入射粒子和靶核看作是點粒子，不考慮其體積。入射時假定每一個入射粒子都以固定的能量垂直射向靶核所在平面?？紤]到靶核質(zhì)量遠大于入射粒子，散射過程中可以認為靶核位置固定。

模擬中，首先應(yīng)進行初始化，設(shè)定散射過程中的幾個參量，主要包括入射粒子能量、入射粒子初始坐標及靶核坐標。入射粒子的坐標與速度確定后，即可對一次散射進行模擬。模擬過程中，認為入射粒子只受庫侖力作用。整個過程可以看作由若干微小時間間隔Δt 組成，每一段間隔內(nèi)，粒子所受庫侖力近似被認為不變。利用特定的算法得出微小時間間隔后入射粒子的速度與坐標。通過多次迭代運算即可得到此次散射過程的軌跡，最后求出該次散射的散射角。通過蒙特卡羅方法隨機給定多個粒子的初始坐標，再利用上述方法模擬出每個粒子的散射軌跡及散射角。對大量粒子的結(jié)果進行統(tǒng)計后即可得到散射幾率。

可以看出在上述模擬流程中，單次散射的軌跡算法對模擬效率起到了關(guān)鍵性作用。因此，下面對多種模擬算法進行詳細的討論和比較。

1.2 軌跡模擬算法

1.2.1 歐拉方法

歐拉方法即為文[8—9]中所使用的算法，是一種以迭代為基本思想的常微分方程的數(shù)值解法。

公式（1）給出了位置函數(shù)r 的Taylor 展開[10]，

其中r(t)是t 時刻α 粒子所在的位置，r(t+Δt)是t 時刻又經(jīng)過一個微小時間段Δt 后α 粒子所在的位置。

取展開式的前兩項，得到：

該公式說明僅需初始時刻的r(t)和t 時刻的r 的微分值（即速度v）即可近似得到t+Δt 時刻的r 值。該算法的誤差為二階。

顯然，歐拉方法的誤差較大，為滿足模擬精度的要求，需要將步長Δt 取的較小，而Δt 過小又會導(dǎo)致模擬效率較低。因此在本研究中嘗試使用其他的算法。

1.2.2 基本Verlet 算法

基本Verlet 算法是一種求解牛頓運動方程的數(shù)值方法，其作為一種積分方法，被廣泛運用于分子動力學(xué)模擬、行星運動等領(lǐng)域。

對位置函數(shù)r 進行Taylor 展開：

可以看出，在該方法中，為了得到t+Δt 時刻的r值，需要t 及t-Δt 時刻的r 值以及t 時刻的加速度。該方法的誤差為四階[11]，精度要高于歐拉算法，但在計算過程中并沒有顯式求出速度，這給最終的散射角計算帶來一定的困難。

1.2.3 Velocity Verlet 算法

在基本Verlet 算法的基礎(chǔ)上，人們又發(fā)展出了Velocity Verlet 算法，該算法也有廣泛的應(yīng)用范圍。

分別對t 和t+2Δt 時刻的位置函數(shù)r 進行Taylor展開并取前3 項：

可以看出該方法的誤差同樣為四階[11-13]，與基本Verlet 算法完全一致，但該算法可以顯式給出每一時刻的速度與位置。

1.2.4 算法比較

為了選擇出最佳算法，本文比較了3 種算法的模擬精度及效率。圖1 中的3 組圖給出了1、5、10 MeV 3 種不同入射粒子能量下3 種算法的模擬誤差。從圖中可以看出：在3 種能量下，歐拉方法的誤差都大于2 種Verlet 算法；隨著入射粒子能量的增加，歐拉方法的誤差變化沒有固定規(guī)律，而2 種Verlet 算法誤差則逐漸減小且算法誤差基本一致。

圖1 不同入射粒子能量下3 種算法誤差比較

表1 給出了3 種算法運行萬次的用時。可以看出：歐拉法與基本 Verlet 算法效率基本相同，Velocity Verlet 算法較另外2 種方法效率略有提高。表2 詳細比較了3 種算法的異同。綜合考慮這些因素可以得到：在軌跡模擬中，Velocity Verlet 算法較其他2 種算法誤差更小，效率更高，也更加方便。因此，選定Velocity Verlet 算法作為改進后的模擬程序的軌跡算法。

表1 3 種算法模擬時間比較

表2 3 種算法綜合比較

1.3 變步長處理

確定好軌跡算法后，模擬效率得到了一定的提高，但還達不到要求。本文希望能進一步提高模擬效率。圖2 給出了α 粒子散射過程中距靶平面不同距離時的受力情況?？梢钥闯?，α 粒子在大部分區(qū)間（距離靶核較遠）內(nèi)受到的庫侖力很小，速度變化很小，因此模擬時步長Δt可以取的較大；而α 粒子距離靶核較近時，庫侖力會急劇增大，該階段對最終的散射軌跡影響很大，因此模擬計算的步長Δt需要取的相對較小，才能保證模擬精度。因此，步長Δt的大小與粒子所受的庫侖力有關(guān)。

為了確定步長 Δt的具體形式，本研究采用了將粒子動量在每個模擬步內(nèi)變化控制在一個固定值這一標準。而粒子動量在每個模擬步內(nèi)變化即為其所受到的沖量。這樣就可以根據(jù)公式（11）求出每一步的步長值?？梢钥闯雒恳徊降牟介L與粒子所受庫侖力F成反比。

圖2 散射過程中F-Z 變化曲線

λ為一固定沖量， 在本研究中其大小取為 6.64×10-23N·s。

對程序進行變步長優(yōu)化后需要對其結(jié)果進行檢驗。而檢驗的標準即為變步長后的模擬精度是否仍在可接受范圍內(nèi)。圖3 是對模擬程序進行變步長處理后，3 個入射粒子能量下散射角誤差隨瞄準距離的變化圖。在3 種入射粒子能量下，散射角絕對誤差全部都不超過1°，可以滿足模擬要求。當(dāng)入射能量較高，瞄準距離較短時，誤差相對較大；其余情況下，誤差都較小且比較穩(wěn)定。

圖3 變步長后不同入射能量下的誤差

2 結(jié)果與分析

經(jīng)過算法篩選、步長動態(tài)處理兩方面的優(yōu)化后，最終給出了基于變步長Velocity Verlet 算法的盧瑟福散射模擬程序。下面將其與原來基于定步長歐拉算法的盧瑟福散射模擬程序進行比較。

圖4 是優(yōu)化后程序與原程序的誤差對比圖。圖4(b)給出了原程序的模擬誤差，可以看出：不同入射能量下的絕對誤差普遍小于0.9°；誤差隨瞄準距離的增大先迅速變大而后減小，在約200 fm 附近出現(xiàn)誤差最大值。圖4(a)給出了優(yōu)化后程序的模擬誤差，可以看出：不同入射能量下散射模擬的絕對誤差同樣都小于0.9°；誤差隨瞄準距離的變化與原程序基本一致，也是先迅速變大而后減小，在200 fm 附近出現(xiàn)最大值。兩張圖對比，可以得到：在較低入射能量下，優(yōu)化程序精度與原程序基本一致；在較高入射能量下，優(yōu)化程序的精度比原程序更好一些。

圖4 優(yōu)化程序、原程序誤差比較

表3 給出了原程序和優(yōu)化后的程序在不同入射能量下運行萬次所需要的時間?？梢钥闯觯谕瑯訔l件下優(yōu)化后的程序運行時間較原程序減小2 個數(shù)量級左右。隨著入射能量的增加，原程序用時與優(yōu)化后程序用時的比值基本保持不變，穩(wěn)定在135 左右。

表3 優(yōu)化前后萬次運行平均用時

本文首先詳細討論比較了多種軌跡算法，并從中優(yōu)選出了Velocity Verlet 算法；然后根據(jù)入射粒子的具體受力情況提出了變步長的模擬方法。通過兩步優(yōu)化可以使得盧瑟福模擬的效率得到大幅提高，優(yōu)化后程序的萬次運行平均用時可以下降到原程序的約1/135。