Программные продукты и системы

Работа посвящена созданию системы моделирования экспериментов по измерению масс заряженных частиц в масс-спектрометре ионного циклотронного резонанса с преобразованием Фурье (Fourier Transform Ion Cyclotron Resonance – FTICR) [1].

FTICR масс-спектрометры обеспечивают высокую точность измерения масс и разрешающую способность, обусловливая их широкое применение в протеомике. Определение масс в FTICR масс-спектрометре основано на резонансном возбуждении частиц в однородном магнитном поле B ионной ловушки на их циклотронной частоте wc=qB/m. Измерив частоту w и зная величину магнитного поля B, можно определить отношение m/q массы m частицы к ее заряду q. Это справедливо лишь для невзаимодействующих частиц, движущихся в отсутствие внешних электрических полей. Как следствие, измеренная в эксперименте частота w отличается от wc, что затрудняет нахождение величины m/q. Одним из способов теоретического исследования влияния не учтенных в формуле параметров на точность определения масс является компьютерное моделирование.

Распространенной программой моделирования в масс-спектрометрическом сообществе является код SIMION (http://www.simion.com), главным образом используемый для нахождения траекторий небольших групп ионов в электрических полях, создаваемых различными конфигурациями электродов.

При реалистическом моделировании FTICR эксперимента необходимо отслеживать эволюцию сотен тысяч взаимодействующих ионов. Чтобы воспроизвести время натурного эксперимента, приходится выполнить миллионы циклов пересчета положений и скоростей частиц. Осуществить подобные вычисления представляется возможным лишь с использованием параллельных компьютеров, что требует создания новых кодов, адаптированных к современным суперЭВМ. Подобный код должен быть интегрирован в привычное экспериментатору окружение.

Целью настоящей статьи является описание системы моделирования FTICR масс-спектромет­ра, построенной на основе разработанного автором на факультете вычислительной математики и кибернетики МГУ имени М.В. Ломоносова (ВМК МГУ) параллельного кода частиц в ячейке PARA­DISE (PARallel ADvanced Ion Simulation Engine).

Математическая постановка задачи

Движение ионов в ловушке масс-спектрометра происходит в соответствии с законами Ньютона и определяется действием силы Лоренца:

здесь mp, qp, ,  – масса p-го иона, его заряд, радиус-вектор и скорость, соответственно. Магнитное поле  предполагается квазистационарным. Действующее на частицу электрическое поле  складывается из удерживающего поля ловушки, радиочастотного возбуждающего поля и кулоновского поля взаимодействия с другими ионами и зарядами, индуцированными на стенках. Оно определяется антиградиентом потенциала: . Если ввести пространственную плотность  заряда, то потенциал  поля внутри ловушки можно искать как решение первой краевой задачи для уравнения Пуассона:

где  Ф/м – электрическая постоянная.

В данной работе рассматривается реальная ловушка c квадратным поперечным сечением. Расстояние между ее торцевыми электродами – a, между парами возбуждающих и детектирующих электродов – d. Постоянный потенциал  прикладывается к торцевым электродам и создает потенциальную яму удерживающего поля. На пару возбуждающих электродов подаются потенциалы . Это радиочастотное электрическое поле возбуждает только те ионы, циклотронная частота wc которых совпадает с частотой . Вследствие резонансного возбуждения частицы начинают двигаться в фазе, образуя ионное облако. Регистрируемым сигналом является разность зарядов, индуцируемых ионными облаками на детектирующих пластинах [1].

Масс-спектром называется Фурье-преобразо­вание этого сигнала – функция, которая показывает вклад в результирующий сигнал ионов с различным значением m/q.

Таким образом, задача состоит в интегрировании уравнений движения ионов и вычислении заряда, индуцированного на детектирующих пластинах ловушки масс-спектрометра.

Метод решения задачи

Вычислительные затраты простейшей схемы «частица–частица» растут квадратично с увеличением числа частиц, поэтому для решения поставленной задачи был выбран метод частиц в ячейке [2]. В этом подходе уравнения движения решаются в области, а для расчета электрического поля вводится сетка. Сеточная плотность заряда находится взвешиванием частиц на сетке. Потенциал определяется из решения уравнения Пуассона. Действующие на частицу силы вычисляются посредством интерполяции по массиву сеточных значений. Детальному описанию алгоритма в применении к рассматриваемой задаче посвящена отдельная публикация [3].

Программная реализация

Так как трехмерный код подразумевает интенсивную работу с многомерными массивами, в основу программной реализации был положен хорошо зарекомендовавший себя в подобных приложениях язык Fortran 90. Чтобы использовать преимущества общей памяти и обеспечить переносимость программы, для распараллеливания были использованы директивы OpenMP (http://www.openmp.org). Изначально программа проектировалась для SMP-системы IBM eServer pSeries 690 Regatta (факультет ВМК МГУ), но с широким распространением многоядерных процессоров появилась возможность запуска кода на настольных рабочих станциях.

Процедуры, требующие распараллеливания, разделяются на три группы: вычисления по частицам, вычисления на сетке, решение уравнения Пуассона.

В циклах вычисления весов интерполяции, распределения заряда на сетку, нахождения действующего на частицу поля и интегрирования уравнений движения итерации являются независимыми по частицам, и для их распараллеливания применена директива omp do. В процессе моделирования некоторые частицы могут навсегда покинуть пределы ловушки. Шаги алгоритма для вылетевших ионов выполнять не нужно, значит вместо статической диспетчеризации omp schedu­le(static) в указанных циклах можно использовать динамическую планировку omp schedule(dynamic, chunk) с определяемой пользователем порцией данных chunk.

При нахождении сеточной плотности заряда имеется опасность возникновения условий гонки, когда в соседних ячейках сетки оказывается более одной частицы и различные нити должны одновременно обновлять значение плотности заряда в общем для этих соседних ячеек узле сетки. В данном случае применение директивы omp atomic допустимо лишь при наличии аппаратной поддержки подобных атомарных операций. Поэтому вместо одного разделяемого нитям выделяются индивидуальные массивы для хранения сеточной плотности. После цикла по частицам проводится суммирование данных из неразделяемых массивов. В случае большого числа узлов имеет смысл распараллелить эту редукцию по третьему измерению и использовать статическую планировку.

Процедура решения уравнения Пуассона использует алгоритм типа Хокни, основанный на комбинации серии независимых двухмерных быстрых преобразований Фурье (БПФ) и метода прогонки. БПФ осуществляется средствами параллельной библиотеки FFTW (www.fftw.org). При распределении циклов метода прогонки применяется статическая диспетчеризация.

Архитектура системы

Система моделирования может использоваться как независимо, так и в связке с натурным экспериментом. Ее архитектура приведена на рисунке. Исходными данными для PARADISE являются состав анализируемого вещества (аналита) и сигнал возбуждения. Количественный зарядово-массовый состав с учетом изотопного распределения генерируется программой emass (www.hel­sinki.fi/science/lipids/software.html), при этом относительные интенсивности изотопных кластеров берутся из экспериментального спектра, а общее число частиц определяется пользователем. Возбуждение осуществляется либо синусоидальным сигналом с линейно меняющейся частотой, либо SWIFT-сигналом. Последний представляет собой обратное Фурье-преобразование от спектральной формы, задаваемой экспериментатором (stored waveform inverse Fourier transform, SWIFT). С помощью графического интерфейса программы AWE экспериментатор выделяет на масс-спектре диапазоны масс (возможно, неперекрывающиеся), которые должны быть подвергнуты возбуждению, и сгенерированный SWIFT-файл подается на вход программы моделирования (при этом AWE используется так же, как и при проведении натурного эксперимента).

Анализ полученного в моделировании ICR-сигнала выполняется программой AWE и разработанной автором утилиты Power. Программа Power использует библиотеку FFTW для вычисления спектра, что дает возможность рассматривать сигналы, число отсчетов в которых не является степенью двойки. Power позволяет анализировать любой выделенный участок сигнала, применять к нему весовые окна, дополнять нулями. В процессе моделирования визуализация осуществляется работающей в системе X Window программой XDraw, а для получения изображений печатного качества подготовлен набор Matlab-скриптов. Для визуализации эволюции ионных облаков предназначена программа IonVis.

Программы AWE и IonVis разработаны сотрудниками Института атомной и молекулярной физики AMOLF (Амстердам, Нидерланды) (www.asms.org/asms06pdf/A062707.pdf).

В заключение отметим, что рассмотренная система используется в Институте энергетических проблем химической физики РАН (г. Москва) и Институте AMOLF для моделирования различных режимов функционирования ионных ловушек масс-спектрометров нового поколения. Система позволяет осуществлять прямое сравнение результатов компьютерного моделирования с натурным экспериментом.

Литература

1.  Marshall A.G., Hendrickson C.L., Jackson G.S. Fourier transform ion cyclotron resonance mass spectrometry: a primer // Mass Spectrom. Rev. 1998. 17, pp. 1–35.

2.  Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987.

3.  Позднеев А.В. Численное моделирование эволюции ионных облаков в масс-спектрометре методом частиц в ячейке // Вестн. Моск. ун-та. 2008. № 3. С. 11–19. (сер. 15. Вычисл. матем. и киберн.).