CC BY
5УДК 541.37; 541.1; 538.975
А. В. Казак, Н. В. Усольцева, В.В. Соцкий
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ПЛАВАЮЩИХ СЛОЕВ АЛКОКСИЗАМЕЩЕННЫХ лезо-ТЕТРАФЕНИЛПОРФИНА И ИХ МЕТАЛЛОКОМПЛЕКСОВ
COMPUTER SIMULATION OF FLOATING LAYERS OF ALKOXYSUBSTITUTED meso-TETRAPHENYLPORPHINE AND THEIR METAL COMPLEXES
Ивановский государственный университет, НИИ наноматериалов 153025 Иваново, ул. Ермака, д. 39. E-mail: alexkazak86@gmail.com
Проведена апробация применения молекулярно-динамического моделирования для многоатомных систем на современных графических контроллерах. Впервые предложена методика численного моделирования плавающих слоев, основанная на применении комбинированных граничных условий для метода молекулярной динамики. По разработанной модели проведены расчеты, согласующиеся с экспериментальными данными для алкоксизамещенных мезо-тетрафенилпорфина.
Ключевые слова: алкоксизамещенные тетрафенилпорфина, моделирование, молекулярная динамика, плавающие слои.
The approbation of the molecular dynamics simulation application for polyatomic systems on modern graphic controllers is made. The method of numerical simulation for floating layers based on the application of combined boundary conditions for the molecular dynamics method is offered for the first time. According to the received model the calculations which were consistent with experimental data for alkoxysubstituted meso-tetraphenylporphines were carried out.
Key words: alkoxysubstituted meso-tetraphenylporphines, computer simulation, molecular dynamics, floating layers.
В настоящее время одной из приоритетных задач науки является получение новых функциональных материалов с заданными свойствами и устройств на их основе. Производные порфиринов, а также их металлокомплексы являются для этих целей одними из наиболее перспективных материалов. Они участвуют в основных процессах обеспечения жизнедеятельности (гемоглобин, миоглобин, цитохром, хлорофилл), нашли применение в качестве катализаторов, полупроводников, материалов для оптических запоминающих устройств, материалов для нелинейной оптики [1]. Что очень важно, эти соединения проявляют высокую химическую и термическую стабильность, а их хорошо известные оптические и эмиссионные свойства позволяют их использовать для научных и практических целей [2].
Наиболее перспективным методом создания устройств на основе данных материалов является метод Ленгмюра — Блоджетт [3]. Основные достоинства данного метода: возможность получения однородных структур молекулярного разрешения, контроль толщины получаемого слоя, создание гибридных наноматериалов путем комбинации слоев различных соединений [3, 4]. Данный метод нашел практическое примене-
© Казак А. В., Усольцева Н. В., Соцкий В. В., 2013.
ние в таких областях техники, как радиоэлектроника и сенсорика. В связи с этим перспективно получение и изучение пленок Ленгмюра — Блоджетт.
Путем введения заместителей в макромолекулу соединения либо введения металла в структуру можно варьировать свойства вещества, что отражается на формировании тонкопленочных материалов. Однако проведение лабораторных экспериментов зачастую связано со значительными временными и материальными издержками. Одним из способов решения данной проблемы может служить создание системы математических моделей, которые позволят получить интересующие характеристики материалов на наноуровне и определить стратегию создания новых молекулярных устройств на их основе. Поэтому в работе с целью апробации применения молекулярно-динамического моделирования плавающих слоев для многоатомных систем на современных графических контроллерах результаты компьютерного эксперимента сравнивались с полученными ранее данными по надмолекулярной упаковке плавающих слоев пятнадцати ал-коксизамещенных мезо-тетрафенилпорфина [5—8].
Материалы и методы исследования
В работе использовались экспериментальные данные, полученные ранее [5—8] для плавающих слоев 15 новых производных мезо-тетрафенилпорфина и их металло-комплексов. Синтез 15 новых соединений с жезо-замещением в молекуле (рис. 1) осуществлялся по методу, описанному ранее [9, 10].
Ленгмюровские слои формировали из растворов соединений I (С = 0,0445 %), II (С = 0,0296 %), III (С = 0,0245 %), IV (С = 0,0249 %), V (С = 0,03597 %), VI (С = 0,0445 %), VII (С = 0,0415 %), VIII (С = 0,0427 %), IX (С = 0,0402 %), X (С = 0,05053 %), XI (С = 0,0481 %), XII (С = 0,043 %), XIII (С = 0,0449 %), XIV (С = 0,0365 %) и XV (С = 0,0441 %) в хлороформе на установке фирмы ФГУП «ГНЦ НИОПИК» (Москва) в широком диапазоне исходных степеней покрытия поверхности с. Скорость сжатия слоя составляла 55 см2/мин [5—8].
\\ ы
I м = 2Н ^ = Н R2 = Н Rз = -ОС4Н9
II м = 2Н ^ = Н R2 = Н Rз = = -ОС16Н33
III м = 2Н R I = —ОС4Н9 R2 = Н Rз = Н
IV м = 2Н = -ОС16Н33 R2 = Н Rз = Н
V м = Си ^ = Н R2 = Н Rз = = -ОС16Н33
VI м = Си = -ОС16Н33 R2 = Н Rз = Н
VII м = 2п ^ = Н R2 = Н Rз = = -ОС16Н33
VIII м = 2п = -ОС16Н33 R2 = Н Rз = Н
IX м = № ^ = Н R2 = Н Rз = = -ОС16Н33
X м = № = -ОС16Н33 R2 = Н Rз = Н
XI м = Pd ^ = Н R2 = Н Rз = = -ОС16Н33
XII м = 2Н ^ = Н R2 = Н Rз = Н
XIII м = Си ^ = Н R2 = Н Rз = Н
XIV м = 2п ^ = Н R2 = Н Rз = Н
XV м = № ^ = Н R2 = Н Rз = Н
Рис. 1. Химическая структура соединений I—XV [5—8]
Математический аппарат
В качестве метода, позволяющего моделировать системы, состоящие из нескольких тысяч частиц, нами был выбран метод молекулярной механики [11]. Основу метода составляет описание исследуемой системы при помощи потенциалов парного взаимодействия. Полная энергия системы складывается из энергии взаимодействия валентных связей, валентных углов, торсионных углов, несвязанных и электростатических взаимодействий:
N
ипот = X кг (г - О2 + К (а - а0)2 +
+ -^(1 + ^(фп-ф0)) +
A B
\
V г12 Г6 ,
+ W 2
S г
Для моделирования системы необходимо иметь минимальный набор констант упругости, равновесные значения длин валентных связей, углов, а также параметры несвязанных взаимодействий и заряды на атомах. Индивидуальный подбор этих параметров для каждой молекулы превращает данный подход в сложную задачу, поэтому в практике моделирования пользуются так называемыми силовыми полями. Они представляют собой набор констант, оптимизированных для определенных классов соединений. Это позволяет с определенной степенью достоверности переносить параметры на исследуемую систему. В зависимости от того, в какую функциональную группу входит атом, ему присваивается определенный тип. По типам атомов, входящих в систему, происходит подстановка параметров силового поля. Также в ходе расчетов необходимо однозначно указать, какой потенциал должен быть выбран для данной пары взаимодействующих частиц. Для системы необходимо составить некое однозначное описание. Мы применили матричное описание, когда для системы составляется матрица связей, в которой указывается какой тип взаимодействия нужно выбирать для данной пары атомов. Заполнение матрицы проводится при помощи графического интерфейса. При этом пользователем определяются начальные координаты и указываются валентно связанные пары атомов.
Для подстановки констант силового поля необходимо провести анализ структуры, при этом:
1) выделяются циклические фрагменты и находится количество атомов, образующих цикл;
2) по ближайшему окружению, типу связи и данным о цикличности подставляются соответствующие типы атомов. Нами использовались типы атомов, соответствующие силовому полю AMBER [12]. На заключительном этапе по типам атомов заполняется матрица силовых постоянных, которая содержит константы упругости, а также равновесные значения при связях, образующих молекулярную систему (рис. 2).
На данном этапе система получает полное механическое описание, но ее структура не соответствует реальной. Для поиска координат атомов применяют методы оптимизации из условия нахождения минимума потенциальной энергии. При построении системы, если неизвестны координаты всех атомов, для нахождения ее структуры используются методы минимизации. Суть их заключается в поиске низкоэнергетических конфигураций молекулярной системы с помощью численных методов поиска экстремумов функций многих переменных. При этом предполагается, что искомая конфигурация находится в области глобального минимума потенциальной энергии [11].
Рис. 2. Анализ молекулярной системы
На практике из-за сложного профиля потенциальной функции поиск глобального минимума превращается в очень сложную вычислительную проблему. Для поиска локальных минимумов потенциальной энергии наиболее широко применяются методы, использующие градиент потенциальной энергии [13, 14]. Мы применили комбинацию градиентных методов оптимизации (наискорейшего спуска и Флетчера-Ривса).
Когда система получает механическое описание и геометрия ее оптимизирована, можно проводить численные эксперименты методом молекулярной динамики. Основу метода составляет численное решение уравнений второго закона Ньютона для системы взаимодействующих частиц:
d2 г (/)
Щ ; ' dt1
или
= Fг (г), г = 1,2,..., N (1)
^ г т—1 = Р.
^ г = 1,2,..., N, (2)
г
где г. - радиус-вектор г -го атома, т. - его масса, Р. - суммарная сила, действующая на -ый атом со стороны остальных частиц, которая определяется как частная производная соответствующего потенциала по координатам, взятая с обратным знаком:
Р (') (3)
ОГ
Здесь: г = {г1, г2,..., гп }; и (г) - потенциальная энергия, зависящая от взаимного расположения всех атомов.
Далее числено решают уравнения движения, вычисляя на каждом шаге все силы, новые скорости и координаты частиц. Для решения этих систем были применены две явные схемы: Эйлера и Рунге — Кутта 2-го порядка [15].
В ходе расчетов важны температурные зависимости. Температура определяется как средняя кинетическая энергия, приходящаяся на одну степень свободы системы:
1 п
Т^) = ^ ^^ . (4)
Если моделирование системы проводится с наложением упругих связей, то N будет равно количеству атомов, из которых состоит система, kB - постоянная Больцмана.
В качестве метода поддержания постоянства температуры системы нами был выбран термостат Берендсена [16]. В данном методе взаимодействие с тепловым резервуаром явно не учитывается, а задается силой трения Qi = т1Л(г )у . Значение коэффициента Л выбирается таким образом, чтобы сила Qi обеспечила изменение энергии системы по закону:
ШЕ 1
= (Е0 - Ек (г)). (5)
ш гЕ
Здесь Е - полная энергия изолированной системы (при отсутствии взаимодействия с резервуаром сохраняется), гЕ - характерное время взаимодействия с резервуа-
N
1, Ек mi (Уi) - кинетическая энергия системы, Е0 = — ^ЬТ0- константа,
ром, Ек = ^ тх yi) - кинетическая энергия системы, е 0 = —
I=1 2
равная средней кинетической энергии, соответствующей температуре резервуара То.
шрг =-?Ц- + Л(г)туг г = 1,2,..., N (6)
ш дх1
Уравнения движения метода имеют вид:
т1Шх1 = р i I = 1,2,..., N (7)
ш
Л =
1
Г Е, ^
2Те
0 -1
(8)
КЕк (О
Данный метод также используется для выбора схемы интегрирования:
АЕтах = Ек-1 - Ек = — (Ео - Ек (г)) г , (9)
ГЕ
здесь АЕтах - максимально возможное изменение энергии системы за промежуток времени г (шаг интегрирования) при использовании термостата. Если изменение энергии в процессе расчета больше этого значения, то считается, что порядок точности интегрирования недостаточен. На основании данного критерия происходит переключение со схемы Эйлера первого порядка точности на схему Рунге — Кутта со вторым порядком точности.
Численные эксперименты, проводимые в неограниченном фазовом пространстве, подходят лишь в том случае, если моделируется система, находящаяся в газовой фазе. Однако моделирование систем, находящихся в жидком или кристаллическом состоянии, требует поддержания в системе постоянных макроскопических характеристик, таких как объем, давление, плотность. Для соблюдения этих условий на систему накладывают ограничения, называемые граничными условиями. Способами ограничения моделируемого фазового пространства являются зеркальные [17] или периодические граничные условия [18].
С применением описанных выше методов нами была разработана система моле-кулярно-динамического моделирования на основе предыдущих разработок [19], позволяющая проводить численное моделирование динамики формирования ленгмюровских пленок.
Реализация расчетов на графических контроллерах
Поскольку метод молекулярной динамики требователен к ресурсам, перевод расчетов на графические контроллеры является перспективным направлением, ввиду высокой производительности данных устройств [20]. Однако графические ускорители являются специализированными устройствами, изначально не предназначенными для решения сложных математических задач. Появление интерфейса программирования СЦОА [21] предложило для вычислений на видеокартах простую и удобную модель разработки программного обеспечения. В то же время особенности архитектуры этих устройств, требуют создания новых алгоритмов [22].
Наиболее затратной операцией при решении задач молекулярной динамики является вычисление суммы сил, действующих на атомы, поскольку количество уравнений пропорционально квадрату количества частиц в системе. Поэтому суммирование лучше проводить на графическом ядре, а интегрирование уравнений движения, занимающих меньшую часть вычислений (количество пропорционально числу частиц), - на центральном процессоре. Функция передачи данных от центрального процессора к графическому контроллеру является синхронизирующим барьером для всех потоков, выполняющихся на видеокарте. Также данная схема наиболее удобна для систем, содержащих несколько графических контроллеров.
Для реализации молекулярной динамики в рамках модели программирования СЦОА была предложена следующая концепция: поскольку для эффективного использования ресурсов графического контроллера требуется создание большого количества потоков, каждому атому выделяется отдельный поток. При этом в нем рассчитываются все действующие на атом силы (рис. 3). Тогда эффективность использования графического ядра растет пропорционально количеству атомов в системе.
Рис. 3. Схема использования графического ядра в задачах молекулярной динамики
Описание молекулярной системы требует большого количества оперативной памяти для хранения координат, параметров атомов, силовых постоянных внутримолекулярного взаимодействия. Большим объемом обладает только глобальная память видеокарты, однако доступ к ней может достигать 400—600 циклов. Частично эта проблема решается использованием схемы расчетов со множеством потоков - запросы на чтение и запись чередуются с расчетами на графическом ядре. В то же время более быстрым доступом обладает константная и общая память. Ограничением для константной памяти является доступ к ней для записи только с CPU, общая же память открыта для записи и чтения графическому ядру. Однако объем общей памяти составляет 16 Кб на один мультипроцессор, а объем данных, используемых в процессе расчета достигает
сотен мегабайт. Для решения этой проблемы был создан алгоритм блочного копирования данных из глобальной в общую память. Каждый поток копирует одно значение из массива глобальной памяти в массив локальной памяти (рис. 4), при этом на N обращений к глобальной памяти приходится N*N-N обращений к общей памяти, за счет чего и обеспечивается прирост производительности программы. Сравнительные численные эксперименты показали, что использование данного алгоритма обеспечивает четырехкратный прирост производительности приложения.
Суммирование взаимодействий по ^ым атомам
Загрузка блока данных в общую память и синхронизация потоков
Рис. 4. Схема работы алгоритма блочного копирования данных
Для сравнения эффективности использованного алгоритма проводились следующие вычислительные эксперименты: запускались приложения скомпилированные с одинаковыми ключами (ключ -02 для оптимизации) и брались молекулярные тестовые системы одинаковой размерности с одинаковыми параметрами моделирования (шаг интегрирования, температура, объем). Затем измерялось количество итераций, выполненное на центральном процессоре и видеокарте за минуту реального времени. Для сравнения использовались процессор Intel Core 2 Duo 2.2 GHz (параллельный вариант программы) и видеокарты Nvidia 9800GT и Nvidia 8600 GT. Зависимость эффективности использования графических ядер от количества частиц представлена на рис. 5.
GeForce 8600GT Intel Core2Duo GeForce 9800GT
О 1000 2000 3000 4000 5000
Число частиц, N
Рис. 5. Зависимость эффективности использования графических ядер от количества частиц в моделируемой системе
Максимальный прирост скорости расчетов в сравнении с двухпроцессорным вариантом (процессор Intel Core 2 Duo 2.2 Ггц) составляет 28,1 раза (видеокарта GeForce 9800 GT) и 49,2 раза (видеокарта GeForce 295 GT).
Результаты исследований
Моделирование плавающих слоев
Моделирование плавающих слоев проводилось в несколько этапов: на первом этапе проводились неэмпирические квантово-химические расчеты методом Хартри — Фока для уточнения парциальных зарядов на атомах моделируемых соединений. Далее полученные заряды подставлялись в уравнения молекулярной динамики.
В качестве модели, выбранной для моделирования воды, использовалась модель TIP3P. Адекватность выходных данных программы и соответствие модели оценивалось по расчету функций радиального распределения g(r) для атомов кислорода. При этом в процессе расчета производилась запись межатомных расстояний, пересчет конечной
dN dN
функции осуществлялся по формуле g (r) =--—, где--количество молекул в
dV 4л r2 р dV
сферическом слое, р - плотность [23]. Моделирование проводилось при температуре 298 К с наложением периодических граничных условий. Конечная функция радиального распределения (рис. 6) соответствует функциям, полученным в работе [24], что подтвердило адекватность получаемых данных в нашей реализации программного обеспечения и возможность его дальнейшего использования на следующих стадиях моделирования.
2,5-,
2,0-
1,5-
1,0 -
0,5-
0,0-
-1-1-'-Г"
2 3 4
~~1-
6
г, А
—1—I—1—I
9 10 11
Рис. 6. Функция радиального распределения O-O для модели TIP3P
Поскольку формирование слоев происходит на поверхности воды, это вносит дополнительные трудности в создание математической модели. Для адекватного моделирования поведения воды необходимо поддерживать макроскопические характеристики, такие как давление и плотность в системе. В то же время для описания поведения молекул вещества на поверхности требуется микроскопический подход (поскольку количество молекул невелико). С целью решения данной проблемы нами впервые был
применен следующий подход: на фазовое пространство, ограничивающее водный слой, накладывались зеркальные граничные условия по осям х, y и z. Для ограничения фазового пространства слоя вещества также применялись зеркальные граничные условия, однако они накладывались только по осям х и z. Таким образом можно поддерживать необходимые макрохарактеристики для объема воды, в то время как перемещение молекул вещества дополнительно не ограничивается по вертикали. Также, за счет раздельных граничных условий можно легко осуществить имитацию сжатия в ленгмю-ровской ванне. Этого можно достигнуть путем изменения верхних граничных условий, не изменяя условий для слоя воды.
Для выявления закономерностей при формировании плавающих слоев проводились численные эксперименты с незамещенным порфирином, тетрафенилпорфином, жезо-алкилоксизамещенными тетрафенилпорфина с заместителями в орто- и пара- положении а также их металлокомплексов (соединения I—XV). Численные эксперименты с незамещенным порфирином проводились для того, чтобы выявить влияние центрального фрагмента на особенности формирования слоев и для изучения влияния заместителей. Изначально молекулы размещались в случайном порядке, с покрытием, соответствующим газовой фазе. Далее размеры граничных условий уменьшались с шагом 1 Ä. Моделирование проводилось до установления постоянных значений по температуре и полной энергии системы (переход системы в равновесное состояние). Типичное время численного эксперимента составляло 500 пс.
Модель слоя после сжатия для соединения III на поверхности воды показана на
рис. 7.
б
а
Рис. 7. Модель слоев после сжатия для соединения III на поверхности воды. Атомы водорода и заместители условно не показаны : а — вид сверху, б — вид сбоку
В ходе расчетов регистрировалась сила, действующая со стороны молекул на стенки расчетной ячейки. По полученным данным были построены зависимости силы от площади покрытия, подобные изотермам п - А сжатия (рис. 8), которые можно качественно сравнивать друг с другом.
0123456789
А, нм2/мол
Рис. 8. Зависимость силы, действующей со стороны молекул на стенки расчетной ячейки, от площади покрытия для соединения III
В ходе численных экспериментов были установлены следующие зависимости: в процесс формирования монослоев существенный вклад вносит взаимодействие центральных фрагментов, в результате которого происходит образование ассоциатов. С увеличением размера молекул (введением заместителей) ассоциация происходит при сжатии слоя. В начальный момент определяющей является сила взаимодействия с поверхностью субфазы. При сжатии (либо увеличении концентрации) расстояния между молекулами уменьшаются, что влечет за собой образование ассоциатов. Введение металлов приводит к увеличению силы взаимодействия между молекулами, в результате чего образуются агрегаты, из которых уже формируется поверхностный слой.
Данные моделирования находятся в хорошем соответствии с результатами ранее проведенного эксперимента [5—8].
Заключение
Разработана система молекулярно-динамического моделирования, позволяющая проводить расчеты на современных программируемых видеокартах.
Впервые предложена методика численного моделирования плавающих слоев, основанная на применении комбинированных граничных условий для метода молекулярной динамики.
По разработанной модели проведены расчеты, результаты которых хорошо согласуются с полученными ранее экспериментальными данными для плавающих слоев пятнадцати алкоксизамещенных мезо-тетрафенилпорфина и их металлокомплексов с замещением (-ОС4Н9 и -ОС16Н33) в пара- или орто-позициях.
Авторы выражают благодарность за предоставленные производные порфина д-ру
хим. наук, проф. Семейкину А. С.
Работа выполнена при финансовой поддержке гранта № 16.740.11.0206 ФЦП Минобр-
науки РФ.
Список использованной литературы
1. Barker C. A., Xianshun Zeng, Bettington S. et al. // Chem. Eur. J. 2007. Vol. 13, № 23. P. 6710 — 6717.
2. Sandeep Kumar // Liq. Cryst. 2009. Vol. 36, № 6/7. P. 607—638.
3. Блинов Л. М. // Успехи физических наук. 1988. Т. 155, № 3. С. 443—480.
4. КовальчукМ. В. // Вестник российской академии наук. 2003. Т. 73, № 5. С. 405.
5. Казак А. В., Усольцева Н. В., Быкова В. В. и др. // Жидкие кристаллы и их практическое использование. 2010. Вып. 4. С. 90—97.
6. Казак А. В., Усольцева Н. В., Юдин С. Г. и др. // Жидкие кристаллы и их практическое использование. 2011. Вып. 2. С. 52—58.
7. Казак А. В., Усольцева Н. В., Юдин С. Г. и др. // Жидкие кристаллы и их практическое использование. 2011. Вып. 3. С. 47—53.
8. Kazak A., Usol'tseva N, Bykova V. et al. // Mol. Cryst. Liq. Cryst. 2011. Vol. 541. P. 28—34.
9. Быкова В. В., Усольцева Н. В., Семейкин А. С. и др. // Жидкие кристаллы и их практическое использование. 2007. Вып. 4. С. 67—74.
10. Быкова В. В., Усольцева Н. В., Семейкин А. С. и др. // Жидкие кристаллы и их практическое использование. 2008. Вып. 2. С. 28—34.
11. Allinger N. L., Burkett U. Molecular Mechanics. ACS, Washington DC. 1982. 326 p.
12. Cornell W. D., CieplakP., Bayly C. I. et al. // J. Am. Chem. Soc. 1995. Vol. 117. P. 5179—5197.
13. Fletcher R. Practical Methods of Optimization. John Wiley & Sons, New York, 1980. 450 p.
14. Gill P. E., Murray W., Wright M. H. Practical Optimization. New York, Academic Press, Inc., 1981.
15. Калиткин Н. Н. Численные методы. М. : Наука, 1979. 512 с.
16. Berendsen H. J. C., Postma J. P. M., van Gunsteren W. F. et al. // J. Chem. Phys. 1984. Vol. 81. P. 3684—3690.
17. Ying-Kui G., Zeng-Yuan G., Xin-GangL. // Chinese Phys. Lett. 2001. Vol. 18. P. 71—73.
18. Балабанян А. Г., Мазо М. А., Ронова И. А. // Химия и компьютерное моделирование. Бутле-ровские сообщения. 2002. № 8. С. 1—12.
19. Соцкий В. В. Свидетельство № 2007613655 об официальной регистрации программы для ЭВМ «Система молекулярно-динамического моделирования LCChem» в Федеральной службе по интеллектуальной собственности, патентам и товарным знакам. М., 2007.
20. Боресков А. В., Харламов А. А. Основы работы с технологией CUDA. Изд-во: ДМК, 2010. 232 с.
21. Документация по интерфейсу программирования CUDA. CUDA Documentation. URL : http://www.nvidia.com/object/cuda_develop.html.
22. Боярченков А. С., Поташников С. И. // Вычислительные методы и программирование. 2009. Т. 10. С. 9—23.
23. Allen M. P., Tildesley D. J. Computer simulation of liquids. New York, Oxford University Press, 1991. 385 p.
24. Mark P., Nilsson L. // J. Phys. Chem. A. 2001. Vol. 105, № 43. P. 9954—9960.
Поступила в редакцию 5.12.2012 г.
