Символьный процессор для расчета импульсных откликов ЯМР двухспиновой системы Текст научной статьи по специальности «Математика»

Вестник Физико-технического института

Крымского федерального университета имени В.И. Вернадского Том 1 (67-69). № 1. 2017. С. 35-49

Journal of Physics and Technology Institute of V.I. Vernadsky Crimean Federal University Volume 1 (67-69). No. 1. 2017. P. 35-49

УДК 537.635; УДК 004.421

СИМВОЛЬНЫЙ ПРОЦЕССОР ДЛЯ РАСЧЕТА ИМПУЛЬСНЫХ ОТКЛИКОВ ЯМР ДВУХСПИНОВОЙ СИСТЕМЫ

Полулях С. Н. *

Физико-технический институт, Крымский федеральный университет имени В. И. Вернадского, Симферополь 295007, Россия *E-mail: sergey. polulyakh@cfuv. ru

Рассматривается реализация символьных вычислений средствами объектно-ориентированного программирования для квантовомеханического расчета намагниченности, формирующейся при импульсном воздействии на систему двух спинов в постоянном магнитном поле с учетом магнитных дипольных взаимодействий. Подробно рассматривается построение объектной модели математического выражения для намагниченности спиновой системы. Приводятся объекты языка программирования для представления спиновых операторов и волновых функций этих операторов, экспоненциальных и тригонометрических функций. Особое внимание уделено программированию операций с дробными коэффициентами, содержащими корень квадратный из двух в знаменателе. В качестве примера приводятся результаты расчета двухимпульсного отклика при произвольных длительностях импульсов, расстройке и величине дипольного взаимодействия. В рамках предложенного подхода удается получить выражение в общем виде для отклика спиновой системы при использовании последовательности, содержащей до четырех импульсов. В частных случаях полученные выражения совпадают с результатами, известными ранее. Ключевые слова: символьные вычисления, магнитный резонанс, двухспиновая система.

PACS: 76.60.-k; PACS: 02.70.-c ВВЕДЕНИЕ

В настоящее время известны специализированные компьютерные программы, обеспечивающие символьные преобразования математических выражений — символьные процессоры. Как правило, в таких программах пользователь вводит исходные математические выражения, придерживаясь правил ввода. Благодаря этому программа различает числа, переменные, функции и т. п. Для математических преобразований выбранного выражения пользователь использует команды меню «разложить на множители», «упростить», «разложить в ряд», и т. п. Символьные процессоры допускают написание пользовательских скриптов, позволяющих программировать сложные вычисления.

Символьные вычисления находят применение в теории магнитного резонанса. Например, в работах [1, 2] сообщается о вычислениях произведения операторов в пакете Mathematica. В работе [3] для символьных вычислений спиновой динамики предложена специализированная программная оболочка SD-CAS, использующая библиотеки компьютерной алгебры YACAS.

В настоящей работе предлагается специализированный символьный процессор для расчета импульсных откликов ЯМР двухспиновой системы, реализованный с помощью методов объектно-ориентированного программирования. Основное внимание в работе уделяется собственно построению объектной модели исходного

математического выражения, используемого для расчета ядерной намагниченности. А именно, выделение в исходном математическом выражении частей, которые могут быть представлены объектами программного кода. Рассматривается программирование правил символьных вычислений с использованием выделенных объектов. Визуализация результатов вычислений осуществляется путем представления объектов в одном из известных текстовых форматов, таких как LaTEX, Html, RTF и т. д.

1. ФИЗИЧЕСКАЯ МОДЕЛЬ

Рассмотрим систему, состоящую из двух магнитных частиц с одинаковым гиромагнитным отношением, находящихся в постоянном магнитном поле (магнитное поле одинаково для каждой из частиц). При этом учитывается магнитное дипольное взаимодействие частиц. Физический пример такой системы — протоны в изолированной молекуле воды [4-6].

Во вращающейся системе координат гамильтониан H двухспиновой системы с учетом только секулярной части дипольных взаимодействий имеет вид (в единицах % = 1, где % — постоянная Планка) [4]

H = -Aö-(/zi + Iz2) + Л-(4• Izi • Iz2-(/+i • I_2 + I_i • I+ 2)). (1)

Здесь Aö — расстройка между частотой вращения системы координат и частотой ядерного магнитного резонанса (ЯМР), определяемой гиромагнитным отношением и величиной постоянного магнитного поля, Л — параметр, характеризующий величину дипольного взаимодействия, Izi, /+г = /хг ± iIyi —

компоненты оператора спина i -й частицы (i = 1,2. Например, Izl — z-компонента оператора спина первой частицы).

При действии на систему импульсов переменного магнитного поля будем полагать, что амплитуда переменного магнитного поля <, выраженная в единицах частоты посредством гиромагнитного отношения, велика как по сравнению с расстройкой < >> A< , так и по сравнению с величиной дипольных взаимодействий < >> Л. Тогда единственное взаимодействие, которое учитывается во время действия импульса — взаимодействие с переменным магнитным полем. Гамильтониан такого взаимодействия (в единицах % = 1) имеет вид [4-6]

R(± x) = +< \/хХ + Iх2 ), (2)

если переменное поле приложено вдоль оси Ox вращающейся системы координат. При направлении поля вдоль положительного направления оси Ox используется

R(+x), R( — при направлении поля вдоль отрицательного направления оси Ox. Если поле приложено вдоль оси Oy, то гамильтониан имеет вид

R±y) = +<1 \ly1 + Iy2 ) . (3)

В настоящей работе рассчитывается макроскопический параметр — поперечная намагниченность M+ = Mx + iMy ансамбля невзаимодействующих друг с другом

двухспиновых систем. Наблюдаемое значение намагниченности вычисляется с помощью оператора матрицы плотности р [4]. Намагниченность в момент времени t есть

M+(t) = Sp(p(t )•(/+1 +1+2 )). (4)

Для оператора матрицы плотности используется уравнение Лиувилля [4]

Р = iP(t),H] = Ш H - H p(t)). (5)

ot

Если в течение времени t гамильтониан H не зависит от времени, то решение уравнения (5) имеет вид1

p(t) = exp(- i • H • t) • p(o) • exp(+ i • H • t). (6)

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

Предполагается, что перед началом действия первого импульса спиновая система находится в состоянии термодинамического равновесия (намагниченность направлена вдоль постоянного магнитного поля). Тогда можно полагать, что в начальный момент времени р = Iz = Izl + Iz2 [4].

Намагниченность ансамбля двухспиновых систем в момент времени t после окончания действия j -го импульса описывается выражением

-iHt „-iRjT, -iHtj-1 -iRj-1T,-1

M +(t) = Sp(e~iHte~iRjTje~iH'J 1 e-^...eI2 x

\ • (7)

х е^..^^-1 е™'-1 )

Здесь тк — длительность к -го импульса, Я — гамильтониан, описывающий взаимодействие спиновой системы с переменным магнитным полем во время действия к -го импульса, tk — длительность временного интервала после к -го импульса (к = 1,2,...,'). Конкретный вид гамильтониана Як определяется выражениями (3) или (4). В дальнейшем временной интервал свободной прецессии после первого импульса обозначается как просто т , а временной интервал после последнего импульса обозначается как t.

1 В этом можно убедиться непосредственной постановкой (6) в (5).

2. ВОЛНОВЫЕ ФУНКЦИИ

Для вычисления следа матрицы в (7) необходимо использовать волновые функции. Для спиновых операторов известны волновые функции, различающиеся магнитным квантовым числом. В случае двухспиновой системы волновую функцию

можно представить в виде | тх , где щ — магнитное квантовое число, относящееся к первой частице. Магнитное квантовое число щ изменяется только спиновыми операторами с индексом 1. Соответственно, щ относится ко второй частице. В случае спина I = 1/2 магнитное квантовое число может принимать значения ± 1/2 . Для краткости записи 1/2 будем опускать, удерживая только знак. Тогда полный набор волновых функций системы двух частиц со спином I = 1/2

каждая представим в виде

1 + +>, 1+"), Ь+>, 1""> .

(8)

Функции (8) являются собственными функциями оператора 12 = + 1г2. Однако, легко видеть, что эти волновые функции не являются собственными для гамильтонианов (1) — (3). С другой стороны, собственные функции этих гамильтонианов можно представить в виде суперпозиции волновых функций (8),

см. Таблицу 1. Собственные функции гамильтониана (1) образуют

ортонормированный базис. Аналогично ортонормированный базис образуют

собственные функции

Р?) и

РР

Используя собственные функции и собственные значения из таблицы 1, выражение (7) представим в виде

(

М + (? ) = £

(Уп1 \е ~'Н1\ У т) {Уп1 \р]п2) {р]П2

Р пМ \Рп2

...X

-V-

ехр (-'Еп )

ехр

(-'42Ъ)

X...

РХп з) \р\з

-'Я1Т1

р1З) рп з

РХп^\Р\ 4

Р\*)\РРп4

ехР ("'4 зг1)

ехР ('4 4Г1)

(9)

X...

Рп5/ \ Рп5

Р]пз)р]п51¥п6)уп6 \е'Н1\¥п6 )&п6 |1+1¥„ 1)

ехр

)

о„

ехР ('Еп6? )

Здесь целые числа п1, п2, пЗ, ... нумеруют волновые функции и соответствующие им собственные значения. Каждое из этих целых чисел принимает

е

с

п\.п2

X

п3п4

е

в

п6,п\

значения от 1 до 4 в соответствии с данными таблицы 1. Суммирование в (9) проводится по всем целым числа п1, п2, п3, ... количество которых зависит от

числа возбуждающих импульсов. В случае импульсной последовательности из N

2 N

импульсов выражение (9) содержит 4 слагаемых. Однако часть этих слагаемых обращается в ноль, а часть слагаемых оказывается подобными для некоторых сочетаний волновых функций.

Таблица 1. Собственные функции и собственные значения гамильтонианов (1) — (3)

Гамильтониан Собственная функция Собственное значение

Н к>=1++> - Аа + Л

^ ^ +">+Т2|'+> - 2Л

0

Аа + Л

Импульс вдоль оси Ох Я{± х} *) = -++>+2+-+--+>+--) -а

=- + +} +- +---- + - ^--) 0

= - ^++)+1|^-++1|-- 0

а

Импульс вдоль оси Оу Я^±у) <) = - ++)+- +->+--+>----> -а

=-+++---+>+--- 0

0

<) = - + + - +-)-- + - ---> а

Верхний индекс у волновых функций Рпу в выражении (9) показывает

порядковый номер возбуждающего импульса, к которому относится волновая функция. В зависимости от направления переменного магнитного поля во время

действия импульса должны использоваться функции

Рп] или

рп ) из таблицы 1.

При переходе от (7) к (9) использована квантово-механическая единица в виде

1=кп>(^п|, 1 =

РХп) РХп

1 =

Рп) [РП

(10)

В выражении (9) можно выделить комплексные числа С и О , получающиеся в результате произведения волновых функций

СП,т = \П \рт) , ОП,т = {Рп \ут) :

(11)

комплексные числа А и В, получающиеся в результате действия операторов /2 = /г1 + /г2 и /+ = /+1 + /+2 на волновые функции

Ап,т = (Р^ ^ Р^) , Вп,ш = \п МО

(12)

и комплексные экспоненты

+1ЕЛ

е

= \п \е±Ш\к >,

±1£1 т

е

1 - ,п}

= Р

е 1

(13)

я,

Р]п

где Еп = \п |Н| у/^ - п -е собственное значение оператора (1), £3п = \ Р,

п -е собственное значение оператора я во время действия -го возбуждающего импульса.

Основная сложность при вычислении намагниченности (9) состоит в большом числе слагаемых. Представляется интересным автоматизировать этот процесс.

3. ОПРЕДЕЛЕНИЕ ИМПУЛЬСНОИ ПОСЛЕДОВАТЕЛЬНОСТИ

В настоящей работе ставится задача получить выражение для временной зависимости намагниченности М+ = М+ ) после последнего возбуждающего импульса согласно выражению (9). Отклик системы - результат действия последовательности импульсов. После каждого импульса следует интервал свободной прецессии. Для задания импульсной последовательности достаточно задать два параметра: порядковый номер импульса в последовательности и направление переменного магнитного поля во время действия импульса. По умолчанию предполагается, что после п -го возбуждающего импульса следует интервал длительностью Тп п+1.

Непосредственно возбуждающий импульс характеризуется углом поворота О , равным произведению длительности импульса т на амплитуду переменного поля

а : О = т • а . В коде программы тип импульса, то есть направление переменного магнитного поля, задается перечислимым типом (рис. 1).

public enum EPulseType {

к = 0, кгл = 1, ^ = 2, Iym = 3

}

Рис. 1. Пример задания типа импульса в синтаксисе C#

Импульсная последовательность представляет собой массив, каждый элемент которого есть тип импульса (перечислимый тип "EPulseType"). Такой объект может быть легко представлен в наглядной форме с помощью строки в RTF, LaTex, Html или другом стандартном формате. Например, RTF представление последовательности из трех импульсов, заданное в тексте программы как EPulseType[] sequence = { EPulseType.Ix, EPulseType.Ixm, EPulseType.Iy}, имеет вид Sequence: a/+X) - x - a2(-X) - T23 - a3(+Y) - t.

4. КОМПЛЕКСНЫЕ ЧИСЛА И МАТРИЧНЫЕ ЭЛЕМЕНТЫ

Компьютерную программу для расчета намагниченности (9) будем строить используя подходы объектно-ориентированного программирования. Достаточно легко определить класс для представления комплексных чисел. Фактически нам достаточно иметь два поля для действительных чисел: действительной и мнимой частей комплексного числа. В этом же классе можно определить операции над комплексными числами: сложение, произведение, модуль числа, проверка на равенство и т. д.

Заметим (см. Таблицу 1), что для СФ гамильтониана (1) числовые множители содержат корень квадратный из двух в знаменателе. При проведении символьных вычислений представляется важным сохранить именно такое представление. Для этого предлагается вместо действительных чисел использовать числа вида

(14)

c + dyl 2

где а, b , c и d — целые числа. С точки зрения программного кода это число может быть представлено объектом с четырьмя целочисленными полями (рис. 2). Простые арифметические операции над такими числами сводятся к сложению и умножению целых чисел и результатом простых арифметических операций над числами вида (14) также является число вида (14). Например, действительное число

представляется объектом с a = 1, Ь = c = 0 и d = 1. Число 1/3 представляется объектом с a = 1, е = 3 и Ь = d = 0 . Преимущество представления (14) состоит в том, что удается удержать точное символьное представление дробных

коэффициентов и коэффициентов, содержащих

42.

public class SNumber {

public int A { get; set; }

public int B { get; set; }

public int C { get; set; }

public int D { get; set; }

/// <summary>

/// Constructor: creates zero /// </summary>

public SNumber() {

A = 0; B = 0; C = 1; D = 0;

}

/// <summary> /// Constructor /// </summary>

public SNumber(int a, int b, int c, int d) {

if (c == 0 && d == 0)

throw new Exception("Error create SNumber"); A = a; B = b; C = c; D = d;

}

/// <summary>

/// Multiplication operator * definition /// </summary>

public static SNumber operator *(SNumber X, SNumber Y) {

int nA = X.A * Y.A + 2 * X.B * Y.B; int nB = X.B * Y.A + X.A * Y.B; int nC = X.C * Y.C + 2 * X.D * Y.D; int nD = X.D * Y.C + X.C * Y.D;

if (nA == 0 && nB == 0) {

nC = 1; nD = 0;

}

return new SNumber(nA, nB, nC, nD);

}

Рис. 2. Пример определения класса для представления чисел в виде (14) и определения операции умножения в синтаксисе C#

Для представления волновых функций достаточно иметь массив из 4-х комплексных чисел. Каждое из чисел есть коэффициент перед собственной функцией оператора Iz . При этом коэффициенты всегда располагаются в строгом порядке: первый элемент массива есть коэффициент перед волновой функцией | + +), второй элемент — коэффициент перед |+—^ и т. д. Операции для нахождения скалярного произведения функций (11) программируются достаточно легко. Например, для скалярного произведения волновых функций можно

определить оператор умножения «*». Если переменная psi представляет волновую функцию , а функция \ представлена переменной А, то скалярное произведение есть комплексное число. Фрагмент С#-кода по определению

класса WaveFn приведен на рисунке 3.

class WaveFn {

/// <summary>

/// Function size

/// </summary>

public const int FN_SIZE = 4;

/// <summary> /// Members /// </summary> public Complex[] C;

/// <summary> /// Default constructor /// </summary>

public WaveFn() {

C = new Complex[FN_SIZE]; for (int pos = 0; pos < FN_SIZE; pos++) C[pos] = new Complex();

}

/// <summary>

/// Multiplication operator * definition /// </summary>

/// <param name="fA">Left wave function</param> /// <param name="fB">Right wave function</param>

public static Complex operator *(WaveFn fA, WaveFn fB) {

Complex c = new Complex();

for (int pos = 0; pos < FN_SIZE; pos++)

c.Add(fA.C[pos].Conjugate * fB.C[pos]); return c;

}

} ..

Рис. 3. Фрагмент С#-кода для определения класса волновой функции

В приведенном примере использован дополнительно определенный класс Complex. Метод «Add» прибавляет к комплексному числу новое комплексное число (аргумент функции). Метод «Conjugate» возвращает комплексно-сопряженное число.

Операторы могут быть представлены перечислимыми типами в силу их конечного числа. Результат действия оператора на волновую функцию есть новая волновая функция. Например,

1г (1 [ - 1 ' ' / ' - 1 ' / ' - 1 ' / ' - 1

2 2 2 2 (15)

1 \ 1 \ 1 \ 1| \ 8 ' 8 '8 ' 81 '

В этом случае исходная волновая функция содержит элементы С[0]=1/2, С[1]= 1/2, С[2]= 1/2, C[3]= 1/2. Функция, полученная в результате действия оператора на исходную функцию, содержит элементы С[0]= 1/8, С[1]= 1/8, С[2]= 1/8 и С[3]= 1/8. Код программы содержит определения для волновых функций, получающихся

в результате действия оператора Iг на волновые функции р¿у и оператора 1+ на

функцию . Этого достаточно для вычисления матричных элементов (13).

Собственно вычисление матричного элемента (13) проводится в два этапа. Вначале оператор действует на правую волновую функцию, а затем находится произведение левой функции на новую функцию, получившуюся в результате действия оператора.

5. ПРЕДСТАВЛЕНИЕ ЭКСПОНЕНЦИАЛЬНЫХ ФУНКЦИЙ

Вначале рассмотрим экспоненциальные функции (13), соответствующие некоторому к -му возбуждающему импульсу. Для каждого множителя вида

рП

e -iRk^k

ркп) = exp(-iekn Tk) (16)

слева от «центрального» оператора 12 в выражении для намагниченности (9) справа от этого оператора найдется множитель

¿Як* к

Ркт

e

pi) = exp(i£km Tk), (17)

где п и т — разные целые числа, Як — гамильтониан, описывающий действие к -

к к

импульса длительностью *к, £п ибт — п -е и т -е собственные значения оператора Як . Таким образом, в коде программы необходимо уметь представлять выражения вида

ехр (- ¿(екп-ект )тк ). (18)

Заметим, что согласно данных таблицы 1 собственные значения типа S^ могут быть представлены амплитудой поля < с соответствующим числовым множителем. Соответственно, разность Skn — S^ в показателе экспоненты (18) также есть произведение числового множителя на <. В свою очередь, произведение длительности импульса на амплитуду есть угол поворота ак = < -Тк. Тогда для представления экспоненты (18) в коде программы достаточно определить комплексный числовой множитель — коэффициент перед схк, тип импульса и его порядковый номер в последовательности.

Рассуждая подобным образом для участков свободной прецессии (в промежутках между возбуждающими импульсами) получим, что для временного интервала тк после k -го импульса получаются экспоненциальные функции вида

exp(-i(En — Em)■ tkk+i )= exp(^k ■ А< • ^+i )■ exp(Yk ■ Л■ h,k+1 ), (19)

где Pk и Yk — комплексные множители. Для представления таких множителей необходимо определить его тип (А< или Л), соответствующий комплексный множитель (Pk или Yk) и порядковый номер в последовательности. Тип экспоненты задается перечислимым типом, комплексный множитель имеет тип Complex, порядковый номер - целое число.

Тогда каждое из слагаемых под знаком суммы в (9) можно представить в виде

Фг = Cr ■ Xr (Tj )■ Xr (А< )■ Xr (лj )■ Xr Tj—1 )■ Xr (а<—1 )■ Xr (Лj—1 )x

X.., Xr (Ti )■ Xr (<■ Xr (Л ).

Здесь под символом r будем понимать полный набор целых чисел n1, n2, n3, ... по которым проводится суммирование, Tj — тип j -го возбуждающего

импульса. В тексте программы слагаемое Фг (20) может быть представлено объектом с двумя полями: множитель (комплексное число Cr) и упорядоченный

массив объектов X , представляющих экспоненциальные множители.

Для экспоненциальных множителей достаточно легко определяется операция сравнения на равенство: два экспоненциальных объекта равны друг другу если они имеют одинаковый тип, равные комплексные множители и одинаковый порядковый номер в последовательности. Операцию сравнения для элементов типа Фг (20) определим без учета множителей C . Два объекта будем полагать равными друг другу, если равны друг другу все экспоненциальные объекты одного типа для каждого из импульсов. Эта операция полезна тем, что при суммировании два

одинаковых объекта типа Фг можно заменить одним с новым числовым множителем Сг.

Для функций типа Xr можно использовать тригонометрическое представление по формуле Эйлера

exp (i а) = cosa + i -sin а. (21)

Тогда представление функций типа Xr в коде программы следует расширить, добавив перечислимый тип функции: экспонента, косинус или синус. Тип функции следует учитывать при проверке на равенство.

6. АЛГОРИТМ ВЫЧИСЛЕНИЙ

Для вычисления отклика согласно (9) будем использовать рекурсию. Начинать рекурсию будем с середины. А именно, перебирая целые числа n и m в интервале от 1 до 4, вначале вычисляем объект типа Ф для первого возбуждающего импульса.

Фо = ( \h\(í)(pÍ () = A,m - exp(-ife ^)

(22)

В качестве входного параметра для соответствующей функции передается новый объект типа Ф с единичным множителем и пустым массивом функций. Результат вычислений — объект, состоящий из комплексного числа А т и функции

ехр (- ¿(е\ - Е1т )*1). Для последующего интервала свободной прецессии, в соответствии с (9), перебирая числа п' и т' в интервале от 1 до 4, находим

Фо = Фо-(Уп\рП)-(рт \¥т)-ехр ^ )- ехр -Л-^ ) (23)

Для этого пересчитываем комплексный множитель

Ап,т Ап,т Сп'п ^тт (24)

и добавляем в массив две новых функции ехр(Д - Аа - ^ 2 ) и ехр(^1 - Л - ^ 2 )

На выходе получаем объект типа Ф и две волновых функции | и \т, |

(точнее, объекты, представляющие эти функции), стоящие «слева» и «справа». Результат вычислений используется как входной параметр для следующего шага рекурсии. Вычисления проводим для всей импульсной последовательности. Результат вычислений для всей последовательности сохраняем как элемент результирующего массива. При добавлении нового элемента в массив проводим проверку на наличие в массиве равного элемента. Если есть равные элементы, то новый элемент не добавляется в результирующий массив, а пересчитывается комплексный множитель имеющегося элемента.

ЗАКЛЮЧЕНИЕ

В качестве примера на рисунке 4 приведены результаты работы программы для последовательности их двух синфазных возбуждающих импульсов. Для визуализации результатов использовано RTF представление. При ах = я/2 и a2 = я такая последовательность обеспечивает формирование эхо Хана. На рисунке 5 приведены результаты для последовательности из двух импульсов, в первом из которых поле направлено вдоль оси Ox, а во втором - вдоль оси Oy вращающейся системы координат. При = я/ 2 и а2 = я/ 2 такая последовательность обеспечивает формирование «солид»-эхо. Для экспоненциальных функций использовано разложение по формуле Эйлера.

Sequence: a/+X) - X - a,2(+X) - t. M+ = Mx + i My:

+2 cos(ai) sin(a2) sin(Ara t) соб(зЛ t)

+2 sin(ai) sin(Ara X) cos(Ara t) ^(зЛ x) cos(зЛ t)

+2 sin(ai) cos(a2) cos(Affl X) sin(Affl t) cos(зЛ X) cos^ t)

-2 sin(ai) cos(a2) cos(Affl X) sin(Affl t) sin(зЛ X) sin(зЛ t)

-2 sin(ai) cos(2a2) sin(Ara X) cos(Ara t) sin^ x) sin(зЛ t)

+2; cos(ai) sin(a2) cos(Ara t) cos(зЛ t)

-2i sin(ai) sin(Ara X) sin(Ara t) cos(зЛ x) cos(зЛ t)

+2i sin(ai) cos(a2) cos(Ara X) cos(Ara t) cos(зЛ x) ^(зЛ t)

-2i sin(ai) cos(a2) cos(Ara X) cos(Ara t) sin(зЛ x) sin(зЛ t)

+2i sin(ai) cos(2a2) sin(Ara X) sin(Affl t) sin^ X) sin(зЛ t)

Mz:

+2 cos(ai) cos(a2)

-2 sin(ai) sin(a2) cos(Ara X) ^(зЛ x)

Рис. 4. Результат расчета намагниченности для последовательности из двух возбуждающих импульсов, приложенных вдоль оси Ox вращающейся системы

координат

Sequence: a1(+X) - Т - a2(+Y) - t.

M+ = Mx + i My:

-2 cos(a1) sin(a2) cos(Ara t) cos(зЛ t)

+2 sin(a1) cos(Ara т) sin(Ara t) cos(зЛ т) cos(зЛ t)

+2 sin(a1) cos(a2) sin(Ara Т) cos(Affl t) ^(зЛ Т) cos(зЛ t)

-2 sin(a1) cos(a2) sin(Ara Т) cos(Affl t) sin(зЛ Т) sin(зЛ t)

-2 sin(a1) cos(2a2) cos(Ara Т) sin(Affl t) sin(зЛ Т) sin(зЛ t)

+2i cos(a1) sin(a2) sin(Ara t) cos(зЛ t)

+2i sin(a1) cos(Ara т) cos(Ara t) cos(зЛ Т) cos(зЛ t)

-2i sin(a1) cos(a2) sin(Ara Т) sin(Ara t) ^(зЛ т) cos(зЛ t)

+2i sin(a1) cos(a2) sin(Affl Т) sin(Affl t) sin(зЛ Т) sin(зЛ t)

-2i sin(a1) cos(2a2) cos(Ara Т) cos(Affl t) sin(3 Л Т) sinß Л t) Mz:

+2 cos(a1) cos(a2)

+2 sin(a1) sin(a2) sin(Ara Т) cos(зЛ т)

Рис. 5. Результат расчета намагниченности для последовательности из двух возбуждающих импульсов, первый из которых приложен вдоль оси Ox, а второй -вдоль оси Oy вращающейся системы координат

Увеличение числа возбуждающих импульсов приводит к быстрому росту числа слагаемых в (9) и, соответственно, к увеличению времени работы программы. Тем не менее, с использованием предложенного алгоритма, наряду с результатами, приведенными на рисунке 4 и рисунке 5, получены отклики для различных трех- и четырехимпульсных последовательностей. Однако эти результаты не приводятся в силу их громоздкости. При аналитических выкладках без использования символьных компьютерных преобразований как правило анализируются случаи относительно простых углов поворота, равных ж или ж/1, отсутствия расстройки

А( или дипольного взаимодействия Л. В случае предложенного символьного процессора результаты в общем виде получаются для произвольных углов поворота намагниченности во время действия возбуждающих импульсов с учетом как расстройки А( , так и дипольного взаимодействия Л. В рамках предложенного подхода несложно реализовать упрощение конечных выражений для конкретных значений используемых параметров.

Список литературы

1. Shriver J. W. NMR Product-Operator Calculations in Mathematica // Journal of Magnetic Resonance. 1991. Vol. 94. Pp. 612-616.

2. Jerschow A. MathNMR: Spin and spatial tensor manipulations in Mathematica // Journal of Magnetic Resonance. 2005. Vol. 176. Pp. 7-14.

3. Filip X., Filip C. SD-CAS: Spin Dynamics by Computer Algebra System // Journal of Magnetic Resonance. 2010. Vol. 207. Pp. 95-113.

4. Сликтер Ч. Основы теории магнитного резонанса. M.: Мир, 1981. 448 c.

5. Дзюба С. А. Основы магнитного резонанса. Новосибирск, 1997, 138 с.

6. Сергеев Н. А., Рябушкин Д. С. Основы квантовой теории ядерного магнитного резонанса. Логос, 2013, 272 с.

SYMBOLIC PROCESSOR TO COMPUTE PULSE NMR RESPONSES OF THE TWO-SPIN SYSTEM Polulyakh S. N. *

Institute of Physics and Technology, V. I. Vernadsky Crimean Federal University, Simferopol 295007, Russia

*E-mail: sergey. polulyakh@cfuv. ru

The implementation of symbolic computations by means of object-oriented programming is reported for the quantum-mechanical calculation of the magnetization, formed by exciting pulses applied to the two-spin system with magnetic dipole interactions in dc magnetic field. The object model of the mathematical expression which gives the spinsystem magnetization is considered in detail. The objects of programming language to represent the spin operators and wave functions of these operators, exponential and trigonometric functions are discussed. Special attention is paid to the manipulations on fractional factors containing the square root of two in the denominator. The results with arbitrary values of pulses duration, detuning and dipole interactions are given for the two-pulse response as an example. Using the proposed approach it is possible to obtain the expressions in a general way for the response of the two-spin system when a sequence include up to four pulses. The obtained expressions match the known expressions for the particular cases.

Keywords: symbolic computations, magnetic resonance, two-spin system.

References

1. J. W. Shriver, Journal ofMagnetic Resonance 94, 612—616 (1991).

2. A. Jerschow, Journal ofMagnetic Resonance 176, 7-14 (2005).

3. X. Filip, C. Filip, Journal ofMagnetic Resonance 207, 95-113 (2010).

4. C. P. Slichter, Principles of magnetic resonance (Springer-Verlag, 1980).

5. S. A. Dzuba, Osnovy magnitnogo resonansa [Fundamentals of magnetic resonance] (Novosibirsk, 1997) [in Russian].

6. N. A. Sergeev, D. S. Ryabushkin, Osnovy kvantovoj teorii magnitnogo resonance [Fundamentals of quantum theory on nuclear magnetic resonance] (Logos, 2013) [in Russian].

Поступила в редакцию 20.02.2017 г. Принята к публикации 06.05.2017 г. Received February 20, 2017. Accepted for publication May 06, 2017