离子阱质谱中离子轨迹算法研究

贺木易，邵睿婷，冯 焱，郭 丹，张虎忠，徐 伟

(1.北京理工大学，北京 100081；2.兰州物理研究所，甘肃 兰州 730000)



离子阱质谱中离子轨迹算法研究

贺木易1，邵睿婷1，冯 焱2，郭 丹1，张虎忠2，徐 伟1

(1.北京理工大学，北京 100081；2.兰州物理研究所，甘肃 兰州 730000)

高精度数值仿真算法是研究离子阱中离子的运动特征、仿真设计与优化离子阱质谱仪器的关键。本研究探讨了2种基于四阶、五阶龙格-库塔算法的数值方法，求解离子在离子阱内的运动方程——Mathieu方程。通过Matlab编程进行离子轨迹仿真实验，计算精度、初值、时间步长等因素对计算结果的影响。结果表明，五阶龙格-库塔算法可以满足离子阱的仿真要求，能实现高精度离子轨迹计算。

离子阱质谱；离子轨迹仿真；Mathieu方程；龙格-库塔算法

质谱作为一种重要的分析方法，已经广泛应用于化学、生物、环境科学、制药行业、空间探测等领域[1]。作为较复杂的大型仪器，质谱仪器的研发周期长、应用成本高，使得仿真技术成了质谱仪器研发、部件性能优化、实验方案设计的理想选择。离子阱质谱仿真技术主要是针对离子阱质量分析器的机械形状与电压施加方式，计算其内部的电场分布，以此研究待测离子(或离子云)在特定的电场分布和气压条件下的运动特性，如离子运动频率、瞬态激发响应等，并根据离子的运动特性表征离子阱的分析性能，如分辨率、准确度、动态范围等[2-7]。

目前报道的离子轨迹仿真系统有ITSIM、ISIS、SIMION等。其中，ITSIM(Ion Trajectory SIMulator)离子轨迹仿真系统由Purdue大学Graham Cooks课题组开发[8]，该系统注重真空流体力学与离子碰撞能量转换过程的仿真，集成了两种离子/分子碰撞模型，对碰撞导致的离子内能转换进行了深入研究。ISIS(Integrated System for Ion Simulation)由Trent大学R.E.March实验室开发[9]，包括一系列离子轨迹计算模块，如MA、FIM、SPQR等。SIMION由Latrobe大学Don C. McGilvery教授研发[10]，是最早的商业化离子轨迹仿真程序，其最初设计主要应用于离子光学系统的设计与仿真，与离子阱设备尚有一定差异。

针对目前离子阱质谱开发设计过程中大量计算资源的要求，一种高精度、具备大规模离子相互作用计算能力的离子阱仿真系统不可或缺，其核心算法是精确求解不同条件下的离子运动轨迹。现有仿真系统的电场与离子轨迹计算精度相对较低，不能捕获离子阱结构设计或加工过程中的微小变形与偏差对其产生的影响，同时长时间离子轨迹仿真的累计误差相对较大，大规模离子反应仿真的运算时间较长，这些都难以满足日益增长的研发需求。

本工作以本课题组前期编写的仿真程序为基础，比较2种基于龙格-库塔算法的数值方法求解Mathieu方程，探讨算法的影响因素，以期获得耗时短、精度高、可靠性好的离子轨迹算法，来实现离子阱中离子轨迹的高精度仿真。

1 实验部分

1.1 离子阱中离子运动方程

离子阱质量分析器具有多种几何形状电极，所形成的电场亦有所差异，但在其中运动的离子具有相似的运动规律。以线性离子阱为例，在xy平面的对称电极上施加的射频(RF)电场为：

φ0=2(U+VcosΩt)

(1)

其中，U为直流分量，V为交流分量，Ω为射频电场角频率。

采用经典力学分析受力，可获得离子在x、y方向上的运动微分方程：

(2)

(3)

其中，m为离子质量，e为离子所带电荷数，r0为离子阱电极距中心的距离。

经简化可得到Mathieu方程形式的运动方程：

(4)

1.2 数值算法

龙格-库塔算法是一种流行的高精度数值方法，最常用的为经典四阶龙格-库塔算法(classicalfourth-orderRKmethod)，其形式为[11]：

(5)

其中：

k1=f(xi，yi)

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

将Mathieu方程化为一阶微分方程组，代入式(5)，得到：

(6)

(7)

其中，u为离子运动位移，v为粒子运动速度。

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

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

消去参数K，得到[12]：

(8)

(9)

每个算式数值均采用经典力学公式积分获得，每一步算法采用经典四阶龙格-库塔系数获得下一步迭代值。

根据龙格-库塔算法的推导规则，可获得高于四阶的方法，且满足不定方程组的系数有无穷多个[13]。一组五阶龙格-库塔系数由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；

消去参数K，得到：

(13)

(14)

同样，每个算式数值均采用经典力学公式积分获得，每一步算法采用五阶龙格-库塔系数获得下一步迭代值。

1.3 仿真实验

仿真程序基于Matlab编写[14]。仿真条件如下：线性离子阱双曲电极距阱中心5 mm，施加射频电场频率为1 MHz，幅度为200 V；仿真离子m/z195，离子初始位置距阱中心4 mm，初始为静止状态；仿真时间2 ms，步长0.05个射频周期。

在x方向，2种算法得到的离子运动轨迹示于图1，线性离子阱中离子运动呈微扰的三角函数周期运动。在无任何操作，且忽略粒子间相互作用力时，离子能量守恒，最大振幅保持在初始位置4 mm，示于图1a。随着仿真时间的增加，算法误差累积，导致离子的最大振幅改变，示于图1b。因此，可以通过计算相对误差值来衡量算法精度。

注：a. 仿真时间0～0.025 ms；b. 仿真时间1.975～2 ms图1 离子轨迹仿真结果Fig.1 The results of ion trajectory simulation

2 结果与讨论

2.1 仿真精度

在离子分子碰撞等仿真过程中，要在很小的空间内准确判断离子的相对位置和速度，就要求更高的算法精度。在短时间内，2种算法误差呈线性增长。分别将2种算法的相对误差随仿真时间的变化拟合曲线，得到四阶误差曲线y=-0.001 66+39.959 37x，R2=0.998 37；五阶误差曲线y=0.000 828 847+2.877 1x，R2=0.832 73，示于图2。仿真进行到2 ms时，四阶算法误差为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

复杂气相离子反应过程通常发生在几十到几百毫秒，短时间的仿真远不能模拟真实的反应进程。随着时间的增长，四阶算法误差快速累积，导致结果超出合理范围。在300 ms仿真过程中，五阶算法离子轨迹包络线示于图3，拟合曲线为y=0.003 997exp(-2.926+x)，R2=0.999 9。由图3可见，五阶算法误差开始时增长较快，然后逐渐变缓，最终趋于平稳。

2.3 初值稳定性

在仿真过程中，通常会根据实际情况设定不同的初值，分别考量初值对2种算法的影响。设离子初始位置为1、2、3、4 mm，仿真时间3 ms，分别拟合相对误差曲线，以曲线斜率作为衡量误差变化的依据，结果示于图4。由图4可见，初值对2种算法的响度误差、稳定性无影响，但随初值的变化，绝对误差随之成比例的变化。

图4 四阶、五阶算法初值稳定性Fig.4 The stability of the initial parameter (position) of 4th and 5th RK algorithm

2.4 步长的选择

数值算法要求选择合理的步长，这对算法稳定性、精度等有着重要的影响。在离子阱仿真中，步长还会影响时域频域的信号处理与分析。减小步长有利于减小计算误差，同时也会导致仿真时间的增长。如上讨论，采用步长为0.05个射频周期时，五阶算法仍不能满足长时间仿真的精度要求，故需要减小步长以提高精度。

五阶算法分别采用0.05、0.02、0.01、0.005个射频周期进行计算，仿真时间3 ms，结果示于图5。随着步长的增大，仿真误差呈指数增加；随着步长的减小，计算时间呈指数增加。因此，在实验中可按仿真需求选择步长。在几百毫秒量级离子轨迹仿真实验中，选择0.01个射频周期步长即可满足仿真精度。

2.5 应用举例

采用优化条件：五阶RK算法，步长为0.01个射频周期，进行100 ms仿真实验。采用双曲电极线性离子阱，电极距离x、y轴5 mm，电极施加频率为1 MHz，幅度为400 V的射频电场；仿真样品为血管紧张素Ⅱ(m/z524，等效半径reff=0.883 nm)。在各种仿真条件下，离子在x方向上的运动轨迹示于图6。

图5 步长对五阶算法仿真精度影响Fig.5 The effect of the step length on 5th RK in ion trajectory simulation

图6a为离子囚禁过程，在不考虑任何离子间相互作用且无操作的条件下，离子在射频电场作用下稳定存储于离子阱中，且保持能量守恒，即离子振幅恒定，仿真误差低于千分之一。图6b为离子检测过程，其操作过程为扫描射频电场幅度，射频初值为400 V，射频扫描速度为10 kV/s。随着射频增大，离子振幅呈被压缩的趋势，当射频达到一定幅度时，离子脱离稳定区，振幅逐渐增大，最终弹射出阱被检测器检测，实现边界激发。图6c为离子冷却过程，是离子与中性气体分子发生碰撞导致能量转移，离子振幅减小，趋近离子阱中心。仿真条件是以1.33×10-3Pa氮气作为缓冲气体，采用hard-sphere中性气体-离子碰撞模型[15]。图6d为离子激发过程，处于离子阱中心的离子通过施加辅助交流电场(AC)获得能量，实现离子碎裂或二次激发等过程。实验施加的AC电场频率为0.316 MHz，扫描速度5 kV/s，在AC作用下，离子能量增高，振幅逐渐增大。以上结果均与经典理论有较高的契合度，从而验证了算法的可靠性。

3 结论

本研究通过数值仿真实验比较了两种基于龙格-库塔算法的数值方法，并将其应用于求解Mathieu方程获得离子阱中离子运动轨迹。实验结果表明，0.01个射频周期步长的5阶龙格-库塔算法具有较高的仿真精度，可在短时间内实现几百毫秒时长的轨迹计算，能满足离子阱中大规模仿真需求。通过仿真实验举例，实现了一些离子阱操作的模拟，验证了算法的可靠性。

[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]朱建共,郭光真. 四极滤质器中离子轨迹的计算方法研究[J]. 厦门大学学报:自然科学版,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

国家自然科学基金(21205005)；科技部仪器专项(2011YQ0900502，2012YQ040140-07)；真空低温技术与物理重点实验室开放基金(ZDK1401)资助

贺木易(1989—)，男(汉族)，北京人，本科，生物医学工程专业。E-mail: fywhhmy@yeah.com

徐 伟(1981—)，男(汉族)，辽宁人，教授，从事分析仪器研究。E-mail: rayxuxu@gmail.com

时间：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