Научная статья на тему 'Возможности информационной технологии Gaussian в моделировании колебательных спектров фосфорорганических соединений (gb-,Gd-, GF-Agents)'

Возможности информационной технологии Gaussian в моделировании колебательных спектров фосфорорганических соединений (gb-,Gd-, GF-Agents) Текст научной статьи по специальности «Физика»

CC BY
174
117
i Надоели баннеры? Вы всегда можете отключить рекламу.

Аннотация научной статьи по физике, автор научной работы — Элькин Михаил Давыдович, Колесникова Ольга Васильевна, Гречухина Оксана Николаевна

На основании модельных неэмпирических расчетов электронной структуры на примере известных фосфорорганических соединений зарина, зомана и циклозарина (GB-,GD-, GF-Agents) показана возможность предсказательных расчетов ИК и КР спектров высокотоксичных фосфорорганических соединений в рамках информационной технологии Gaussian.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

The article describes the possibility of predictable calculations based on non-empiric calculations of electronic structure on the example of well-known organophosphorous compounds sarin, soman, cyclosarin (GB-, GD-, GF-Agents).

Текст научной работы на тему «Возможности информационной технологии Gaussian в моделировании колебательных спектров фосфорорганических соединений (gb-,Gd-, GF-Agents)»

УДК 539.194

М.Д. Элькин, О.В. Колесникова, О.Н. Гречухина

ВОЗМОЖНОСТИ ИНФОРМАЦИОННОЙ ТЕХНОЛОГИИ GAUSSIAN В МОДЕЛИРОВАНИИ КОЛЕБАТЕЛЬНЫХ СПЕКТРОВ ФОСФОРОРГАНИЧЕСКИХ СОЕДИНЕНИЙ (GB-,GD-, GF- AGENTS)

На основании модельных неэмпирических расчетов электронной структуры на примере известных фосфорорганических соединений -зарина, зомана и циклозарина (GB-,GD-, GF-Agents) показана возможность предсказательных расчетов ИК и КР спектров высокотоксичных фосфорорганических соединений в рамках информационной технологии Gaussian.

M.D. Elkin, O.V. Kolesnikova, O.N. Grechuhina

TECHNIQUE GAUSSIAN AND MODELING OF VIBRATIONAL SPECTRA

OF GB-, GD-, GF-, AGENTS

The article describes the possibility of predictable calculations based on non-empiric calculations of electronic structure on the example of well-known organophosphorous compounds - sarin, soman, cyclosarin (GB-, GD-, GF-Agents).

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

Первый подход, заявленный в работе [1], опирается на известный фрагментарный метод [2] и использует библиотеку силовых постоянных изученных молекулярных фрагментов. До недавнего времени это был доминирующий теоретический метод ИК спектроскопии. На то были веские причины, главная из которых - недостаточная точность имеющихся квантовых методов в оценке таких параметров адиабатического потенциала, определяющего физические и химические свойства соединения. Уязвимым местом подхода являлась методика сшивки фрагментов, неизбежный произвол при формировании базы данных фрагментарных и силовых постоянных и электрооптических параметров. Подавляющее большинство указанных молекулярных параметров было получено путем решения обратных колебательных задач [3]. Для этого требовалось наличие полного набора экспериментальных данных по частотам и интенсивностям колебаний изотопозамещенных соединений при условии, что интерпретация спектров проведена надежно. Возникающие на этом пути трудности хорошо известны и подробно описаны в диссертации [3]. Тем не менее, для ряда соединений и фрагментов органической химии удалось получить оценку гармонических силовых постоянных и электрооптических параметров, что может служить ориентиром при построении структурно-динамических моделей вновь синтезируемых соединений.

Второй подход связан с привлечением неэмпирических квантовых методов расчета структуры и спектров многоатомных молекул. В настоящее время этот подход становится

доминирующим в предсказательных расчетах колебательных спектров сложных молекулярных соединений, а получаемые результаты, представленные, к примеру, в работах [3-4], весьма обнадеживающие. Особенно это касается квантовых методов, связанных с использованием DFT (методов функционала плотности) подхода [5], реализованного в виде известной компьютерной технологии «GAUSSIAN».

Критические замечания в его адрес, высказанные в публикации [1], справедливы и учитываются разработчиками при создании очередной версии программного обеспечения. Основанием этому служит динамика развития данного направления в теоретической колебательной спектроскопии. В программных комплексах последнего поколения появилась возможность анализа ангармонизма колебаний [5].

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

В данной работе, на примере известных высокотоксичных фосфорорганических соединений - зомана, зарина и циклозарина, имеющих общий токсичный фрагмент P-(F,O, CH3), - показана возможность применения неэмпирических квантовых методов для предсказательных расчетов геометрической структуры и колебательных спектров фосфорорганических соединений с точностью, достаточной для спектральной идентификации этих соединений. Расчеты спектров осуществлены в ангармоническом приближении теории молекулярных колебаний. Предпочтение отдано неэмпирическому квантовому методу DFT/B3LYP / 6-31(6-311)G*(**) информационной технологии Gaussian [5]. Это позволит выяснить влияние корреляционных и поляризационных эффектов квантовой химии на поведение спектральных полос.

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

Поскольку мы не располагаем достаточными литературными данными по структуре и спектрам исследуемых фосфорорганических соединений, то естественным является подход, связанный с моделированием адиабатического потенциала соединений, исходя из возможных конформационных моделей и применяя различные базисные наборы в методе функционала плотности [5].

Для молекулярных фрагментов, составляющих фосфорорганические соединения, мы ограничились четырьмя базисными наборами 6-31G*(**) и 6-311G*(**). Два первых набора хорошо зарекомендовали себя при расчете парафиновых, непредельных и ароматических углеводородов [4], остальные - в пробных расчетах ряда фосфорорганических соединений [7].

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

Для атомов углерода и фосфора имеет место гибридизация БР3. Атом кислорода образует мостик между фрагментами парафиновых углеводородов и фенольного кольца с фрагментом, центральный атом которого фосфор. Для парафиновых углеводородов и монозамещенных бензола геометрические параметры и частоты фундаментальных колебаний известны [6]. Они хорошо воспроизводятся неэмпирическими квантовыми расчетами в рамках метода функционала плотности, если осуществлять расчеты в ангармоническом приближении теории молекулярных колебаний. Это подтверждается расчетами колебательных спектров табуна, зомана, зарина и бензофенона соответственно, представленными в работах [7, 8]. Однако указанные расчеты следует считать предварительным этапом построения структурно-динамических моделей

рассматриваемых соединений. Анализ влияния базиса квантового метода, внутреннего вращения отдельных фрагментов и их взаимного расположения на результаты теоретического исследования колебательных спектров в указанных работах не проводился. Такой анализ необходим, поскольку экспериментальными данными в низкочастотной области спектра (ниже 400 см-1) мы не располагаем.

Модельный гамильтониан, выбранный для решения ангармонической колебательной задачи во втором порядке теории возмущения, имеет следующий вид [6]

И"1 = 1 ]+РМ‘0ь +[0а0ъ0с ]+ ]}, (1)

где Таъ - контравариантный метрический тензор; О - естественные колебательные координаты; Ра - соответствующие им операторы импульсов; ГаЪ - квадратичные; ¥аЪс -кубические; ¥аъсй - квартичные силовые постоянные.

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

Е '1 = ш, (у + я, /2) + X ,, (V, +1/ 2) (у, +1/ 2)(1 +1/28, (2)

где ш, (см-1) - частоты гармонических колебаний; х,г (см-1) - поправки ангармонического приближения; V, - квантовые числа колебательного состояния. Везде предполагается суммирование по индексам правой части соотношений.

Явный вид выражений для коэффициентов х„, являющихся функциями гармонических частот колебаний, ангармонических силовых и кинематических постоянных, резонансных выражений вида 1/(ш, ±ш, ±шг) приведен в монографии [4].

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

Согласно проведенным расчетам выбор базиса и модели конформера несущественно влияют на длины связей общего фрагмента зарина, зомана и циклозарина: СО=1.46 - 1.47А, ОР=1.6 - 1.62А, РО=1.47 - 1.48А, РБ=1.60 -1.61А, РС=1.8 -1.81А. Разброс в значении угла мостика СОР составляет 120.3 - 126.9° и определяется моделью

конформера. Изменения остальных углов зависят от базиса и лежат в диапазонах:ОРО=114.1 - 118°, ОРБ=101.8 -102.8°, 0РС=100.9 - 106.7°; БРО=11.6 - 114.3°, С0Р=116.2 - 118.6°, БРС=100.9 - 102.7°. Наибольший разброс имеет место для зарина. Отметим, что в исходных моделях все валентные углы для атомов с гибридизацией БР3 полагались тетраэдрическими.

Расчетные значения двухгранных углов между плоскостью мостика СОР и плоскостями тетраэдрического фрагмента, содержащего атом фосфора, определяются как исходной моделью конформера, так и выбранным базисом для различных молекул. При задании исходных данных одна из связей РХ (Х=Б, О, С) полагалась компланарной плоскости мостика СОР, а остальные связи лежали в плоскостях, полученных поворотом на 120° относительно исходной. Исходная компланарность нарушается более чем на 9°. Изменение углов между плоскостями ОРХ и ОРУ (Х, У= Б, О, С) от исходного значения ё 120° достигает 15° лишь в зарине для одной из моделей в базисе 6-311О. В остальных молекулах остается в пределах 5°, а для метильных групп не превосходит 1,5°. Отступление от тетраэдрических значений валентных углов метильных групп лежит в пределах 1°.

Интересными представляются результаты моделирования геометрической структуры циклозарина. Исходные три модели определяются положением мостика СОР относительно плоскости фенильного кольца. Углы между плоскостями п, п/2, п/4 соответственно. С учетом описанных выше конформаций фрагмента, содержащего центральный атом фосфора относительно связи СО, имеем 32 исходные модели.

Оптимизация геометрии приводит к трем конформерам: одному для случая, когда в исходной модели компланарны связи СО, ОР, РО(РБ), двум - в случае компланарности связей СО, ОР, РС. На длинах валентных связей и валентных углов обоих фрагментов это практически не сказывается (!^~0.005А; Др~1.2°) и является результатом накопления ошибок в численных методах оптимизации геометрии в неэмпирических расчетах методами функционала плотности. Для первой конформационной модели угол между плоскостью мостика СОР и фенильным кольцом лежит в интервале 41.3°-42.6°. Для двух других конформеров: 48.1°-49.1° (Модели п, базис 6-31 О) и 67.5°-67.9° (Модели п/2, п/4, базис 6-31 О): 44.0°-44.7° (Модели п, базис 6-311 О) и 66.1°-67.8° (Модели п/2, п/4, базис 6-311 О).

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

Рассчитанные минимумы адиабатических потенциалов для конформеров зарина, зомана, циклозарина и табуна попадают соответственно в интервалы: -750.18 —750.33; -828.81 - -828.98; -863.29 - -863.46; -798.48 - -798.64 а.е. К понижению значения минимума адиабатического потенциала приводят переход к базису 6-311О и учет поляризационных эффектов.

Частоты фундаментальных колебаний второго фрагмента, отвечающего за токсичные свойства соединений, следует считать результатами предсказательного расчета в табл. 1-3, поскольку приведенные нами экспериментальные данные [7] являются весьма ограниченными и представлены в диапазоне 600-4000 см-1 в большинстве случаев лишь спектрограммами. Характер поведения интенсивности полос в ИК спектрах для одинаковых молекулярных фрагментов парафиновых углеводородов указывает на характеристичность соответствующих колебаний для всех рассматриваемых соединений и хорошо согласуется с экспериментом. Специфичным является и характер спектра второго фрагмента. Здесь легко идентифицируются валентные колебания связей РО, РК, КС, РН. Что касается деформационных колебаний этого фрагмента, то их можно считать характеристическими по частотам и формам колебаний, однако использование этих полос для идентификации соединений проблематично ввиду слабой интенсивности в спектрах ИК и КР.

Анализ внутреннего вращения отдельных фрагментов относительно мостика РОС показывает на несущественное смещение фундаментальных состояний. Это смещение лежит в пределах ~ 10-15 см-1 и, согласно оценкам из работы [1], на процедуре идентификации соединений по их колебательным спектрам не сказывается. Однако интенсивность ряда полос парафиновых фрагментов изменяется существенно. Именно по ним можно судить о конформации метильных групп относительно мостика РОС. Более подробный анализ внутреннего вращения в соединениях данного класса является предметом дальнейшего исследования авторов.

Выводы

1. Неэмпирические расчеты колебательных спектров приведенных фосфорорганических соединений указывают на возможность их привлечения для предсказательных расчетов конформационной структуры и молекулярных параметров в колебательных спектрах и других молекул, принадлежащих к классу фосфорорганических соединений. При этом предпочтение следует отдать методу 6-311О(**). Квантовые расчеты позволяют оценить систему гармонических силовых постоянных, определяющих форму адиабатического потенциала.

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

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

4. Конформация фрагмента, содержащего атом фтора, относительно парафинового фрагмента заметно сказывается на расчетных значениях интенсивностей в низкочастотном диапазоне (ниже 600 см-1).

5. Изменение длин валентных связей и значений валентных углов для парафинового фрагмента находится в границах, приведенных в монографии [6] для парафиновых углеводородов, и не сказывается на положении полос валентных и деформационных колебаний алкильных групп.

6. Идентификация исследуемых соединений может опираться на полосы, отнесенные к колебаниям связей Р-О, Р=О, Р-С, Р-Б. Интенсивность соответствующих полос значительна (табл. 1). Идентификация конкретного соединения производится по полосам, присущим алкильным группам и фенильному кольцу.

Таблица 1

Интерпретация колебательного спектра фрагмента COP-CH3FO в молекуле зомана (С7H12FO2P)

Форма колебаний Vexp V™ ИК КР

Qp=o 1276 1272 - 1282 168 - 197 4.16 - 4.82

Qco 1016 991 - 1007 244 - 308 3.71 - 4.47

Qop - 960 - 971 315 - 365 1.41 - 4.54

Qpf 840 796 - 832 82.5 - 114 1.92 - 3.92

Qpc - 713 - 724 25.3 - 37.9 12.3 - 17.5

Popo - 447 - 489 17.5 - 29.7 2.49 - 4.61

Popc - 425 - 447 11.6 - 37.2 0.56 - 2.44

aopf - 386 - 402 9.51 - 14.6 1.99 - 2.93

- 248 - 275 0.15 - 0.74 0.09 - 0.6

popf - 228 - 234 0.31 - 1.03 0.18 - 0.97

pcop - 124 - 140 6.5 - 8.79 0.30 - 0.57

Примечание. Частоты колебаний (^) в см-1; Иитенсивности в спектрах ИК в Км/моль; интенсивности в спектрах КР А4/аеш; обозначения форм колебаний соответствуют таковым, принятым в работе [6]

Таблица 2

Интерпретация фундаментальных колебаний фрагмента ООР-ОИзРО в циклозарине (Оптимизированная модель 1)

Форма колебаний 6-31 О* 6-31 О** 6-311 О* 6-311 О**

vm ИК КР Vm ИК КР Vm ИК КР Vm ИК КР

Ор=о 1274 142 5.54 1272 135 5.39 1270 154 6.28 1268 145 6.23

Оор 941 168 1.02 938 342 1.14 936 338 931 406 1.29

ОрГ 836 90.6 2.30 834 86.7 2.33 804 128 2.05 803 127.4 1.97

Орс 751 37.3 2.53 751 37.9 2.55 745 36.9 2.20 744 40.3 2.16

рорГ 446 41.8 4.67 446 42.4 4.72 445 47.7 4.98 445 47.7 4.98

роро 412 32.6 3.13 412 30.4 2.89 413 19.2 1.40 413 15.6 1.11

Рфо 371 7.90 1.79 370 7.70 1.76 366 8.55 2.07 366 8.50 2.03

Рссо 309 3.76 1.37 309 3.85 1.38 309 3.97 1.35 309 3.96 1.34

рсро 258 3.19 2.78 259 3.29 2.79 258 3.56 2.84 258 3.66 2.85

Рорс 222 0.93 0.96 222 0.93 0.99 225 1.00 1.31 225 0.99 1.31

Рсор 111 5.99 0.77 112 6.10 0.73 108 6.61 0.71 108 6.66 0.70

Таблица 3

Интерпретация фундаментальных колебаний фрагмента ООР-ОИ3РО в циклозарине (Оптимизированные модели 2 и 3)

Форма колебаний 6-31 О*(**) 6-311О*(**)

Модель п2, п/4 Модель п Модель п2, п/4 Модель п

Vm ИК КР Vm ИК КР Vm ИК КР Vm ИК КР

Ор=о 1292 120 4.68 1278 128 5.33 1293 281 9.73 1274 148 6.57

Оор 936 302 1.84 942 395 2.18 937 229. 1.66 938 383 1.52

ч— р О 819 48,9 3.03 846 82.4 3.78 792 123 1.75 819 63.3 0.89

Орс 746 2,63 1.82 750 12.4 1.75 742 10.4 0.89 751 37.3 2.53

рорГ 455 30,1 1.52 465 53.4 3.99 457 31.9 1.54 464 60.7 4.40

о р ч— 411 26,4 2.04 414 29.1 2.13 409 30.6 2.27 412 34.5 2.38

Рссо 336 9,93 1.28 325 6.36 1.30 335 11.6 1.25 324 5.34 1.23

Роро 317 5,88 4.07 303 5.52 1.75 319 5.25 4.09 298 7.70 2.14

Рорс 257 1,26 1.51 264 2.58 2.17 259 1.91 1.62 262 2.75 1.96

Рсро 242 0,53 1.24 240 0.06 1.50 243 0.61 1.52 243 0.05 1.74

Рсор 102 0,22 1.92 114 2.98 1.44 172 0.10 0.28 178 0.30 0.79

Таблица 4

Предсказательный спектр фундаментальных колебаний зарина (С4Н10РО2Р)

Форма ^ехр Vh Vanh ИК КР Форма ^"ехр Vh Vanh ИК КР

СО I О оа 1487 1538 1481 5.8 3.3 Ррсн - 940 908 31 4.2

СО I О оа 1468 1525 1470 5.1 25 Оее 867 896 871 8.7 7.4

со I О оа 1468 1516 1460 0.1 16 °рр - 847 821 96 2.1

со I О оа 1460 1511 1459 1.2 12 °РС - 770 769 18 2.7

со I О оа 1452 1492 1438 4.1 13 о е О - 723 684 21 21

Ренэ 1452 1491 1434 7.2 16 Роро - 499 475 31 5.8

Ренэ 1390 1449 1399 10 3.9 весе - 475 461 16 0.6

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

РеН3 1372 1437 1381 15 2.9 Рорр 424 429 8.1 1.7

РеН3 1335 1399 1357 3.9 4.5 весе 375 404 395 15 1.7

Роен - 1395 1336 11 11 аорс - 370 354 12 0.8

Ppch - 1382 1328 35 1.7 Росс - 307 302 3.3 1.1

Qpo 1292 1306 1312 179 4.2 xx - 260 264 1.4 0.4

Pccc 1152 1215 1181 10 2.2 aFPc - 253 258 0.7 0.7

Росе - 1172 1134 12 3.7 Popc - 239 236 2.3 1.1

Рссн - 1145 1108 46 2.7 xx 207 222 222 0.1 0.1

Qop - 1010 989 500 3.4 xx 173 146 0.1 0.1

Ppch 974 963 936 78 0.9 Pcop 142 141 8.8 0.5

Pcch 923 960 931 4.4 1.5 xx 64 56 1.4 0.1

Ppch - 950 926 1.3 2.2 xx 33 33 1.6 0.2

* Примечание. Приведены частоты фундаментальных колебаний соответствующих фрагментов из монографии [6]

ЛИТЕРАТУРА

1. Мясоедов Б.Ф. Фрагментарные методы расчета ИК спектров фосфорорганических соединений / Б.Ф. Мясоедов, Л. А. Грибов, А.И. Павлючко // Журнал структурной химии. 2006. Т. 47. № 1. С. 449-456.

2. Грибов Л. А. Методы и алгоритмы вычислений в теории колебательных спектров молекул / Л.А. Грибов, В.А. Дементьев. М.: Наука, 1981. 356 с.

3. Березин К.В. Квантово-механические модели и решение на их основе прямых и обратных спектральных задач для многоатомных молекул: дис. ... доктора физ.-мат. наук / К.В. Березин. Саратов, 2004. 246 с.

4. Пулин В.Ф. Исследование динамики молекулярных соединений различных классов / В.Ф. Пулин, М.Д. Элькин, В.И. Березин. Саратов: СГТУ, 2002. 546 с.

5. Caussian 03, Revision B.03 / M.J. Frisch, G.W. Trucks, H.B. Schlegel et al. Pittsburg PA, 2003. 680 р.

6. Свердлов Л.М. Колебательные спектры многоатомных молекул / Л.М. Свердлов, М. А. Ковнер, Е.П. Крайнов. М.: Наука, 1970. 550 с.

7. Элькин П.М. Методы оптической физики в экологическом мониторинге фосфорорганических соединений / П.М. Элькин, В.Ф. Пулин, А.С. Кладиева // Вестник Саратовского государственного технического университета. 2007. № 2 (25). С. 177-182.

8. Элькин П.М. Электронная структура и колебательные спектры бензофенона / П.М. Элькин, В.Ф. Пулин, И.И. Гордеев // Вестник Саратовского государственного технического университета. 2007. № 2 (25). С. 51-56.

Элькин Михаил Давыдович -

доктор физико-математических наук, профессор кафедры «Информатика»

Саратовского государственного технического университета

Колесникова Ольга Васильевна -

ассистент кафедры «Информатика»

Саратовского государственного технического университета

Гречухина Оксана Николаевна -

ассистент кафедры «Прикладная математика и информатика»

Астраханского государственного технического университета

Статья поступила в редакцию 07.09.07, принята к опубликованию 15.01.08

i Надоели баннеры? Вы всегда можете отключить рекламу.