賀木易，邵睿婷，馮 焱，郭 丹，張虎忠，徐 偉

(1.北京理工大學(xué)，北京 100081；2.蘭州物理研究所，甘肅 蘭州 730000)

?

離子阱質(zhì)譜中離子軌跡算法研究

賀木易1，邵睿婷1，馮 焱2，郭 丹1，張虎忠2，徐 偉1

(1.北京理工大學(xué)，北京 100081；2.蘭州物理研究所，甘肅 蘭州 730000)

高精度數(shù)值仿真算法是研究離子阱中離子的運(yùn)動(dòng)特征、仿真設(shè)計(jì)與優(yōu)化離子阱質(zhì)譜儀器的關(guān)鍵。本研究探討了2種基于四階、五階龍格-庫(kù)塔算法的數(shù)值方法，求解離子在離子阱內(nèi)的運(yùn)動(dòng)方程——Mathieu方程。通過Matlab編程進(jìn)行離子軌跡仿真實(shí)驗(yàn)，計(jì)算精度、初值、時(shí)間步長(zhǎng)等因素對(duì)計(jì)算結(jié)果的影響。結(jié)果表明，五階龍格-庫(kù)塔算法可以滿足離子阱的仿真要求，能實(shí)現(xiàn)高精度離子軌跡計(jì)算。

離子阱質(zhì)譜；離子軌跡仿真；Mathieu方程；龍格-庫(kù)塔算法

質(zhì)譜作為一種重要的分析方法，已經(jīng)廣泛應(yīng)用于化學(xué)、生物、環(huán)境科學(xué)、制藥行業(yè)、空間探測(cè)等領(lǐng)域[1]。作為較復(fù)雜的大型儀器，質(zhì)譜儀器的研發(fā)周期長(zhǎng)、應(yīng)用成本高，使得仿真技術(shù)成了質(zhì)譜儀器研發(fā)、部件性能優(yōu)化、實(shí)驗(yàn)方案設(shè)計(jì)的理想選擇。離子阱質(zhì)譜仿真技術(shù)主要是針對(duì)離子阱質(zhì)量分析器的機(jī)械形狀與電壓施加方式，計(jì)算其內(nèi)部的電場(chǎng)分布，以此研究待測(cè)離子(或離子云)在特定的電場(chǎng)分布和氣壓條件下的運(yùn)動(dòng)特性，如離子運(yùn)動(dòng)頻率、瞬態(tài)激發(fā)響應(yīng)等，并根據(jù)離子的運(yùn)動(dòng)特性表征離子阱的分析性能，如分辨率、準(zhǔn)確度、動(dòng)態(tài)范圍等[2-7]。

目前報(bào)道的離子軌跡仿真系統(tǒng)有ITSIM、ISIS、SIMION等。其中，ITSIM(Ion Trajectory SIMulator)離子軌跡仿真系統(tǒng)由Purdue大學(xué)Graham Cooks課題組開發(fā)[8]，該系統(tǒng)注重真空流體力學(xué)與離子碰撞能量轉(zhuǎn)換過程的仿真，集成了兩種離子/分子碰撞模型，對(duì)碰撞導(dǎo)致的離子內(nèi)能轉(zhuǎn)換進(jìn)行了深入研究。ISIS(Integrated System for Ion Simulation)由Trent大學(xué)R.E.March實(shí)驗(yàn)室開發(fā)[9]，包括一系列離子軌跡計(jì)算模塊，如MA、FIM、SPQR等。SIMION由Latrobe大學(xué)Don C. McGilvery教授研發(fā)[10]，是最早的商業(yè)化離子軌跡仿真程序，其最初設(shè)計(jì)主要應(yīng)用于離子光學(xué)系統(tǒng)的設(shè)計(jì)與仿真，與離子阱設(shè)備尚有一定差異。

針對(duì)目前離子阱質(zhì)譜開發(fā)設(shè)計(jì)過程中大量計(jì)算資源的要求，一種高精度、具備大規(guī)模離子相互作用計(jì)算能力的離子阱仿真系統(tǒng)不可或缺，其核心算法是精確求解不同條件下的離子運(yùn)動(dòng)軌跡。現(xiàn)有仿真系統(tǒng)的電場(chǎng)與離子軌跡計(jì)算精度相對(duì)較低，不能捕獲離子阱結(jié)構(gòu)設(shè)計(jì)或加工過程中的微小變形與偏差對(duì)其產(chǎn)生的影響，同時(shí)長(zhǎng)時(shí)間離子軌跡仿真的累計(jì)誤差相對(duì)較大，大規(guī)模離子反應(yīng)仿真的運(yùn)算時(shí)間較長(zhǎng)，這些都難以滿足日益增長(zhǎng)的研發(fā)需求。

本工作以本課題組前期編寫的仿真程序?yàn)榛A(chǔ)，比較2種基于龍格-庫(kù)塔算法的數(shù)值方法求解Mathieu方程，探討算法的影響因素，以期獲得耗時(shí)短、精度高、可靠性好的離子軌跡算法，來(lái)實(shí)現(xiàn)離子阱中離子軌跡的高精度仿真。

1 實(shí)驗(yàn)部分

1.1 離子阱中離子運(yùn)動(dòng)方程

離子阱質(zhì)量分析器具有多種幾何形狀電極，所形成的電場(chǎng)亦有所差異，但在其中運(yùn)動(dòng)的離子具有相似的運(yùn)動(dòng)規(guī)律。以線性離子阱為例，在xy平面的對(duì)稱電極上施加的射頻(RF)電場(chǎng)為：

φ0=2(U+VcosΩt)

(1)

其中，U為直流分量，V為交流分量，Ω為射頻電場(chǎng)角頻率。

采用經(jīng)典力學(xué)分析受力，可獲得離子在x、y方向上的運(yùn)動(dòng)微分方程：

(2)

(3)

其中，m為離子質(zhì)量，e為離子所帶電荷數(shù)，r0為離子阱電極距中心的距離。

經(jīng)簡(jiǎn)化可得到Mathieu方程形式的運(yùn)動(dòng)方程：

(4)

1.2 數(shù)值算法

龍格-庫(kù)塔算法是一種流行的高精度數(shù)值方法，最常用的為經(jīng)典四階龍格-庫(kù)塔算法(classicalfourth-orderRKmethod)，其形式為[11]：

(5)

其中：

k1=f(xi，yi)

k4=f(xi+h，yi+k3h)

將Mathieu方程化為一階微分方程組，代入式(5)，得到：

(6)

(7)

其中，u為離子運(yùn)動(dòng)位移，v為粒子運(yùn)動(dòng)速度。

K1=vn，L1=f(Tn，un，vn)，f=-(a+2qcosAT)u；

K4=vn+hL3，L4=f(Tn+h，un+hK3，vn+hL3);

消去參數(shù)K，得到[12]：

(8)

(9)

每個(gè)算式數(shù)值均采用經(jīng)典力學(xué)公式積分獲得，每一步算法采用經(jīng)典四階龍格-庫(kù)塔系數(shù)獲得下一步迭代值。

根據(jù)龍格-庫(kù)塔算法的推導(dǎo)規(guī)則，可獲得高于四階的方法，且滿足不定方程組的系數(shù)有無(wú)窮多個(gè)[13]。一組五階龍格-庫(kù)塔系數(shù)由Butcher(Butcher fifth-order RK method)給出，算法形式為[11]：

(10)

其中：

k1=f(xi,yi)

將Mathieu方程化為一階微分方程組，代入式(10)，得到：

(11)

(12)

K1=vn，L1=f(Tn，un，vn)，f=-(a+2qcosAT)u；

消去參數(shù)K，得到：

(13)

(14)

同樣，每個(gè)算式數(shù)值均采用經(jīng)典力學(xué)公式積分獲得，每一步算法采用五階龍格-庫(kù)塔系數(shù)獲得下一步迭代值。

1.3 仿真實(shí)驗(yàn)

仿真程序基于Matlab編寫[14]。仿真條件如下：線性離子阱雙曲電極距阱中心5 mm，施加射頻電場(chǎng)頻率為1 MHz，幅度為200 V；仿真離子m/z195，離子初始位置距阱中心4 mm，初始為靜止?fàn)顟B(tài)；仿真時(shí)間2 ms，步長(zhǎng)0.05個(gè)射頻周期。

在x方向，2種算法得到的離子運(yùn)動(dòng)軌跡示于圖1，線性離子阱中離子運(yùn)動(dòng)呈微擾的三角函數(shù)周期運(yùn)動(dòng)。在無(wú)任何操作，且忽略粒子間相互作用力時(shí)，離子能量守恒，最大振幅保持在初始位置4 mm，示于圖1a。隨著仿真時(shí)間的增加，算法誤差累積，導(dǎo)致離子的最大振幅改變，示于圖1b。因此，可以通過計(jì)算相對(duì)誤差值來(lái)衡量算法精度。

注：a. 仿真時(shí)間0～0.025 ms；b. 仿真時(shí)間1.975～2 ms圖1 離子軌跡仿真結(jié)果Fig.1 The results of ion trajectory simulation

2 結(jié)果與討論

2.1 仿真精度

在離子分子碰撞等仿真過程中，要在很小的空間內(nèi)準(zhǔn)確判斷離子的相對(duì)位置和速度，就要求更高的算法精度。在短時(shí)間內(nèi)，2種算法誤差呈線性增長(zhǎng)。分別將2種算法的相對(duì)誤差隨仿真時(shí)間的變化擬合曲線，得到四階誤差曲線y=-0.001 66+39.959 37x，R2=0.998 37；五階誤差曲線y=0.000 828 847+2.877 1x，R2=0.832 73，示于圖2。仿真進(jìn)行到2 ms時(shí)，四階算法誤差為7.83%，五階算法誤差為0.66%，由此可見，五階算法精度比四階算法高約11倍。

圖2 算法誤差曲線Fig.2 The relative error of 4th and 5th RK algorithm in ion trajectory simulation

2.2 收斂性

圖3 五階算法收斂性曲線Fig.3 The convergence of 5th RK algorithm in ion trajectory simulation

復(fù)雜氣相離子反應(yīng)過程通常發(fā)生在幾十到幾百毫秒，短時(shí)間的仿真遠(yuǎn)不能模擬真實(shí)的反應(yīng)進(jìn)程。隨著時(shí)間的增長(zhǎng)，四階算法誤差快速累積，導(dǎo)致結(jié)果超出合理范圍。在300 ms仿真過程中，五階算法離子軌跡包絡(luò)線示于圖3，擬合曲線為y=0.003 997exp(-2.926+x)，R2=0.999 9。由圖3可見，五階算法誤差開始時(shí)增長(zhǎng)較快，然后逐漸變緩，最終趨于平穩(wěn)。

2.3 初值穩(wěn)定性

在仿真過程中，通常會(huì)根據(jù)實(shí)際情況設(shè)定不同的初值，分別考量初值對(duì)2種算法的影響。設(shè)離子初始位置為1、2、3、4 mm，仿真時(shí)間3 ms，分別擬合相對(duì)誤差曲線，以曲線斜率作為衡量誤差變化的依據(jù)，結(jié)果示于圖4。由圖4可見，初值對(duì)2種算法的響度誤差、穩(wěn)定性無(wú)影響，但隨初值的變化，絕對(duì)誤差隨之成比例的變化。

圖4 四階、五階算法初值穩(wěn)定性Fig.4 The stability of the initial parameter (position) of 4th and 5th RK algorithm

2.4 步長(zhǎng)的選擇

數(shù)值算法要求選擇合理的步長(zhǎng)，這對(duì)算法穩(wěn)定性、精度等有著重要的影響。在離子阱仿真中，步長(zhǎng)還會(huì)影響時(shí)域頻域的信號(hào)處理與分析。減小步長(zhǎng)有利于減小計(jì)算誤差，同時(shí)也會(huì)導(dǎo)致仿真時(shí)間的增長(zhǎng)。如上討論，采用步長(zhǎng)為0.05個(gè)射頻周期時(shí)，五階算法仍不能滿足長(zhǎng)時(shí)間仿真的精度要求，故需要減小步長(zhǎng)以提高精度。

五階算法分別采用0.05、0.02、0.01、0.005個(gè)射頻周期進(jìn)行計(jì)算，仿真時(shí)間3 ms，結(jié)果示于圖5。隨著步長(zhǎng)的增大，仿真誤差呈指數(shù)增加；隨著步長(zhǎng)的減小，計(jì)算時(shí)間呈指數(shù)增加。因此，在實(shí)驗(yàn)中可按仿真需求選擇步長(zhǎng)。在幾百毫秒量級(jí)離子軌跡仿真實(shí)驗(yàn)中，選擇0.01個(gè)射頻周期步長(zhǎng)即可滿足仿真精度。

2.5 應(yīng)用舉例

采用優(yōu)化條件：五階RK算法，步長(zhǎng)為0.01個(gè)射頻周期，進(jìn)行100 ms仿真實(shí)驗(yàn)。采用雙曲電極線性離子阱，電極距離x、y軸5 mm，電極施加頻率為1 MHz，幅度為400 V的射頻電場(chǎng)；仿真樣品為血管緊張素Ⅱ(m/z524，等效半徑reff=0.883 nm)。在各種仿真條件下，離子在x方向上的運(yùn)動(dòng)軌跡示于圖6。

圖5 步長(zhǎng)對(duì)五階算法仿真精度影響Fig.5 The effect of the step length on 5th RK in ion trajectory simulation

圖6a為離子囚禁過程，在不考慮任何離子間相互作用且無(wú)操作的條件下，離子在射頻電場(chǎng)作用下穩(wěn)定存儲(chǔ)于離子阱中，且保持能量守恒，即離子振幅恒定，仿真誤差低于千分之一。圖6b為離子檢測(cè)過程，其操作過程為掃描射頻電場(chǎng)幅度，射頻初值為400 V，射頻掃描速度為10 kV/s。隨著射頻增大，離子振幅呈被壓縮的趨勢(shì)，當(dāng)射頻達(dá)到一定幅度時(shí)，離子脫離穩(wěn)定區(qū)，振幅逐漸增大，最終彈射出阱被檢測(cè)器檢測(cè)，實(shí)現(xiàn)邊界激發(fā)。圖6c為離子冷卻過程，是離子與中性氣體分子發(fā)生碰撞導(dǎo)致能量轉(zhuǎn)移，離子振幅減小，趨近離子阱中心。仿真條件是以1.33×10-3Pa氮?dú)庾鳛榫彌_氣體，采用hard-sphere中性氣體-離子碰撞模型[15]。圖6d為離子激發(fā)過程，處于離子阱中心的離子通過施加輔助交流電場(chǎng)(AC)獲得能量，實(shí)現(xiàn)離子碎裂或二次激發(fā)等過程。實(shí)驗(yàn)施加的AC電場(chǎng)頻率為0.316 MHz，掃描速度5 kV/s，在AC作用下，離子能量增高，振幅逐漸增大。以上結(jié)果均與經(jīng)典理論有較高的契合度，從而驗(yàn)證了算法的可靠性。

3 結(jié)論

本研究通過數(shù)值仿真實(shí)驗(yàn)比較了兩種基于龍格-庫(kù)塔算法的數(shù)值方法，并將其應(yīng)用于求解Mathieu方程獲得離子阱中離子運(yùn)動(dòng)軌跡。實(shí)驗(yàn)結(jié)果表明，0.01個(gè)射頻周期步長(zhǎng)的5階龍格-庫(kù)塔算法具有較高的仿真精度，可在短時(shí)間內(nèi)實(shí)現(xiàn)幾百毫秒時(shí)長(zhǎng)的軌跡計(jì)算，能滿足離子阱中大規(guī)模仿真需求。通過仿真實(shí)驗(yàn)舉例，實(shí)現(xiàn)了一些離子阱操作的模擬，驗(yàn)證了算法的可靠性。

[1] MARCH R E, TODD J F. Quadrupole ion trap mass spectrometry[M]. John Wiley & Sons, 2005.

[2] PLASS W R, LI H, COOKS R G. Theory, simulation and measurement of chemical mass shifts in RF quadrupole ion traps[J]. International Journal of Mass Spectrometry, 2003, 228(2): 237-267.

[3] HAN S J, SHIN S K. Space-charge effects on Fourier transform ion cyclotron resonance signals: Experimental observations and three-dimensional trajectory simulations[J]. Journal of the American Society for Mass Spectrometry, 1997, 8(4): 319-326.

[4] PARKS J H, SZ?KE A. Simulation of collisional relaxation of trapped ion clouds in the presence of space charge fields[J]. The Journal of Chemical Physics, 1995, 103(4): 1 422-1 439.

[5] NIKOLAEV E N, HEEREN R, POPOV A M, et al. Realistic modeling of ion cloud motion in a Fourier transform ion cyclotron resonance cell by use of a particle-in-cell approach[J]. Rapid Communications in Mass Spectrometry, 2007, 21(22): 3 527-3 546.

[6] LONDRY F A, ALFRED R L, MARCH R E. Computer simulation of single-ion trajectories in Paul-type ion traps[J]. Journal of the American Society for Mass Spectrometry, 1993, 4(9): 687-705.

[7] DING L, SUDAKOV M, KUMASHIRO S. A simulation study of the digital ion trap mass spectrometer[J]. International Journal of Mass Spectrometry, 2002, 221(2): 117-138.

[8] BUI H A, GRAHAM COOKS R. Windows version of the ion trap simulation program ITSIM: A powerful heuristic and predictive tool in ion trap mass spectrometry[J]. Journal of Mass Spectrometry, 1998, 33(4): 297-304.

[9] FORBES M W, SHARIFI M, CROLEY T, et al. Simulation of ion trajectories in a quadrupole ion trap: A comparison of three simulation programs[J]. Journal of Mass Spectrometry, 1999, 34(12): 1 219-1 239.

[10]LOCK C M, DYER E W. Simulation of ion trajectories through a high pressure radio frequency only quadrupole collision cell by SIMION 6.0[J]. Rapid Communications in Mass Spectrometry, 1999, 13(5): 422-431.

[11]CHAPRA S C, CANALE R P. Numerical methods for engineers[M]. New York: McGraw-Hill Higher Education, 2010: 727-737.

[12]朱建共,郭光真. 四極濾質(zhì)器中離子軌跡的計(jì)算方法研究[J]. 廈門大學(xué)學(xué)報(bào):自然科學(xué)版,2005,44(1):41-43.

ZHUJian gong, GUO Guangzhen. Studies of calculation methods about ions trace in the quadrupole mass filter[J]. Journal of Xiamen University (Natural Science), 2005, 44(1): 41-43(in Chinese).

[13]BUTCHER J C, MOIR N. Experiments with a new fifth order method[J]. Numerical Algorithms, 2003, 33(1/2/3/4): 137-151.

[14]XIONG X, XU W, FANG X, et al. Accelerated Simulation study of space charge effects in quadrupole ion traps using GPU techniques[J]. Journal of the American Society for Mass Spectrometry, 2012, 23(10): 1 799-1 807.

[15]GUAN S, LI G Z, MARSHALL A G. Effect of ion-neutral collision mechanism on the trapped-ion equation of motion: A new mass spectral line shape for high-mass trapped ions[J]. International Journal of Mass Spectrometry and Ion Processes, 1997, 167: 185-193.

Study of the Simulation Method about Ion Trajectory

HE Mu-yi1, SHAO Rui-ting1, FENG Yan2, GUO Dan1, ZHANG Hu-zhong2, XU Wei1

(1.BeijingInstituteofTechnology,Beijing100081,China;2.LanzhouInstituteofPhysics,Lanzhou730000,China)

A high precision numerical simulation algorithm is the key of simulation design, which is a powerful tool to study the characteristics of the movement of ions in ion trap and the optimization of ion trap mass spectrometer. Two numerical methods based on fourth-order and fifth-order Runge-Kutta algorithm to solve Mathieu equation were discussed. By using the Matrix Laboratory (Matlab) programming, the research studied the impacts of different parameters, such as calculation accuracy, initial value of position and time step, on the calculation results of ion trajectories. The results show that the fifth-order Runge-Kutta algorithm can meet the requirements of simulation and achieve high precision ion trajectory simulation.

ion trap mass spectrometry; ion trajectory simulation; Mathieu equation; Runge-Kutta algorithm

2014-05-22;

2014-06-25

國(guó)家自然科學(xué)基金(21205005)；科技部?jī)x器專項(xiàng)(2011YQ0900502，2012YQ040140-07)；真空低溫技術(shù)與物理重點(diǎn)實(shí)驗(yàn)室開放基金(ZDK1401)資助

賀木易(1989—)，男(漢族)，北京人，本科，生物醫(yī)學(xué)工程專業(yè)。E-mail: fywhhmy@yeah.com

徐 偉(1981—)，男(漢族)，遼寧人，教授，從事分析儀器研究。E-mail: rayxuxu@gmail.com

時(shí)間：2014-12-02；

http:∥www.cnki.net/kcms/doi/10.7538/zpxb.youxian.2014.0070.html

O657.63

A

1004-2997(2015)03-0217-06

10.7538/zpxb.youxian.2014.0070