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

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

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

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

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

1 模擬方法

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

1.1 模擬過(guò)程

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

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

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

1.2 軌跡模擬算法

1.2.1 歐拉方法

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

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

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

取展開式的前兩項(xiàng)，得到：

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

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

1.2.2 基本Verlet 算法

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

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

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

1.2.3 Velocity Verlet 算法

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

分別對(duì)t 和t+2Δt 時(shí)刻的位置函數(shù)r 進(jìn)行Taylor展開并取前3 項(xiàng)：

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

1.2.4 算法比較

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

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

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

表1 3 種算法模擬時(shí)間比較

表2 3 種算法綜合比較

1.3 變步長(zhǎng)處理

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

為了確定步長(zhǎng) Δt的具體形式，本研究采用了將粒子動(dòng)量在每個(gè)模擬步內(nèi)變化控制在一個(gè)固定值這一標(biāo)準(zhǔn)。而粒子動(dòng)量在每個(gè)模擬步內(nèi)變化即為其所受到的沖量。這樣就可以根據(jù)公式（11）求出每一步的步長(zhǎng)值。可以看出每一步的步長(zhǎng)與粒子所受庫(kù)侖力F成反比。

圖2 散射過(guò)程中F-Z 變化曲線

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

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

圖3 變步長(zhǎng)后不同入射能量下的誤差

2 結(jié)果與分析

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

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

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

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

表3 優(yōu)化前后萬(wàn)次運(yùn)行平均用時(shí)

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