Таблица 5.2 Заряды на атомах в молекулах фенола и бензонитрила

Таблица 5.2

Заряды на атомах в молекулах фенола и бензонитрила

(по Лёвдину – для неэмпирических расчетов)

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

Сначала необходимо рассчитать геометрию молекулы, выбрав предварительно базис. Затем в том же базисе вычисляется гессиан, для чего в группе команд $CONTRL указывается параметр RUNTYP=HESSIAN.

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

5.6. Пример. Расчет ИК-спектра
молекулы формальдегида

Цели расчета: Определение критериев выбора базиса и теоретического уровня для расчета колебательного спектра молекулы. Знакомство с результатами расчета вторых производных (гессиана).

Рассчитаем геометрию молекулы формальдегида в нескольких базисах, и в тех же базисах – гессиан. В последнем случае в конце выходного файла содержатся частоты колебаний, а также термохимические данные, которые вычисляются в приближении идеального газа с учетом нулевых колебаний и приводятся к температуре 298.15 K (MP2/6-311G(d,p)):

Основные типы колебаний молекулы формальдегида (рис. 5.3) можно наглядно посмотреть с помощью программы MacMolPlt (меню Subwindow – Frequencies и меню View – Animate Mode).

Расчет дает 12 частот колебаний молекулы, из которых 3 имеют нулевые интенсивности, т.е. не активны в ИК-спектре. Реальный спектр обычно снимают в диапазоне частот 4000-400 см^-1, и мы рассмотрим только частоты колебаний, попадающие в эту область и имеющие отличные от нуля интенсивности (в случае формальдегида таких частот 6, табл. 5.3).

	Валентные асимметричные колебания связей С-Н	νas (C-H)
	Валентные симметричные колебания связей С-Н	νs (C-H)
	Валентные колебания связи С=О	ν (C=О)
	Деформационные колебания угла НСН (с участием связи С=О)	δ (HCH)
	Деформационные колебания угла НСO	δ (HCO)

Р и с. 5.3. Типы колебаний молекулы формальдегида

5.7. Работа №3. Расчет барьеров вращения

Для расчета барьера вращения проводят ряд расчетов с фиксированным торсионным углом. В результате определяется приблизительная геометрия и энергия переходного состояния, соответствующая максимуму на энергетической кривой. Точно переходное состояние можно рассчитать (для сложных структур) путем задания параметра RUNTYP=SADPOINT и включения в исходный файл матрицы вторых производных в группу команд $HESS, которая располагается после группы $DATA. Для этого предварительно нужно рассчитать гессиан с приблизительной геометрией переходного состояния. Расчет гессиана переходного состояния и исходных реагентов дает вклады колебательных и вращательных поправок к электронной энергии частиц при температуре 298.15К. Эти поправки позволяют вычислить энтальпию и энтропию активации, а также константу скорости мономолекулярной реакции.

5.8. Пример. Вычисление барьера вращения

1) построить график зависимости энергии молекулы от диэдрального угла С-С=С-С;

2) вычислить термохимические параметры E-, Z-2-бутена и переходного состояния;

3) вычислить энергию активации и константу скорости изомеризации.

Исходный файл для расчета Е-2-бутена:

$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE NZVAR=30 MAXIT=30 MULT=1 COORD=ZMT $END

$SYSTEM TIMLIM=525600 MEMORY=10000000 $END

$BASIS GBASIS=N31 NGAUSS=6 $END

$GUESS GUESS=HUCKEL $END

$SCF DIRSCF=.TRUE. $END

$STATPT OPTTOL=0.0001 NSTEP=40 IFREEZ(1)=6,9,12 $END

$DATA

Title

C

C 1 1.5055654

C 3 1.4999549 2 122.5719520 1 180.0000000 0

H 2 1.0773805 3 123.4646696 1 180.0000000 0

H 3 1.0859564 2 122.5599595 4 180.0000000 0

H 4 1.0819998 3 114.2174549 5 -42.4221205 0

H 4 1.0877600 7 108.7171911 3 132.7547045 0

H 4 1.0904689 8 106.8201524 7 -113.8758803 0

H 1 1.0815736 2 108.9778867 5 -10.1163130 0

H 1 1.0793932 10 107.9179480 2 114.7055970 0

H 1 1.0878453 10 109.6943955 11 118.3573448 0

$END
В этом файле присутствует параметр IFREEZ(1)=6,9,12 в группе команд $STATPT, который "замораживает" 6-ю, 9-ю и 12-ю координаты (выделены жирным). Для "замораживания" необходимо также в группе команд $CONTRL указать общее количество координат в параметре NZVAR=3M-6 – для нелинейных молекул или NZVAR=3M-5 – для нелинейных молекул. Проведя несколько оптимизаций геометрии при разных фиксированных значениях торсионного угла, строят график зависимости энергии от данного угла. В результате получают следующую зависимость (табл. 5.4, рис. 5.4).
Разности энергий:

ΔE_ZE = 0.0024049341 а.е. = 6.31 кДж/моль;

ΔE_Z≠ = 0.1635117504 а.е. = 429.30 кДж/моль;

ΔE_E≠ = 0.1659166845 а.е. = 435.61 кДж/моль.

Z-изомер E-изомер

Переходное состояние соответствует величине торсионного угла С-С=С-С 90^о, и в данном случае можно не проводить уточненный расчет переходного состояния. Вычисление гессиана (RUNTYP=HESSIAN) для изомеров и переходного состояния дает соответствующие термодинамические параметры колебательных и вращательных состояний частиц при температуре 298.15 К:

KJ/MOL KJ/MOL KJ/MOL J/MOL-K J/MOL-K J/MOL-K

ELEC. 0.000 0.000 0.000 0.000 0.000 0.000

TRANS. 3.718 6.197 -41.197 12.472 20.786 158.962

ROT. 3.718 3.718 -26.670 12.472 12.472 101.924

VIB. 274.369 274.369 256.926 40.815 40.815 58.503

TOTAL 281.806 284.284 189.059 65.758 74.073 319.389
Е-изомер E H G CV CP S

KJ/MOL KJ/MOL KJ/MOL J/MOL-K J/MOL-K J/MOL-K

ELEC. 0.000 0.000 0.000 0.000 0.000 0.000

TRANS. 3.718 6.197 -41.197 12.472 20.786 158.962

ROT. 3.718 3.718 -26.287 12.472 12.472 100.638

VIB. 273.791 273.791 259.612 40.727 40.727 47.555

TOTAL 281.228 283.706 192.128 65.670 73.984 307.154
Переходное состояние

E H G CV CP S

KJ/MOL KJ/MOL KJ/MOL J/MOL-K J/MOL-K J/MOL-K

ELEC. 0.000 0.000 0.000 0.000 0.000 0.000

TRANS. 3.718 6.197 -41.197 12.472 20.786 158.962

ROT. 3.718 3.718 -26.717 12.472 12.472 102.079

VIB. 265.948 265.948 260.036 29.240 29.240 19.828

TOTAL 273.385 275.864 192.123 54.183 62.498 280.869

ΔH^≠(398.15K) = ΔE_(≠-E) + ΔH_(≠-E) = (435.61 - 8.42) = 427.19 кДж/моль;

ΔS^≠(398.15K) = -38.5 Дж/(моль ^. K).
Для превращения Z-изомера в Е-изомер получают энтропию и энтальпию активации:
ΔH^≠(398.15K) = ΔE_(≠-Z) + ΔH_(≠-Z) = (429.30 – 7.84) = 421.46 кДж/моль;

ΔS^≠(398.15K) = -26.3 Дж/(моль ^. K).

где k_B = 1.3806 ^. 10^-23 Дж/K – постоянная Больцмана,

(k_B/h = 2.08358 ^. 10¹⁰ K^-1c^-1).

Возьмем среднее значение для активационных параметров, так как их разность невелика: ΔH^≠ = 424.3 кДж/моль, ΔS^≠ = - 32 Дж/(моль ^. K). Получается следующее значение константы скорости (реакция практически не идет):

.

Насколько близка эта величина к экспериментальной, если экспериментальное значение барьера вращения для алкенов лежит в интервале 250-270 кДж/моль? В чем состоит основная ошибка расчета? Как более точно вычислить барьер вращения? Сравните энтальпийный и энтропийный вклады в барьер вращения. Насколько необходимо вычислять гессиан и соответствующие поправки при проведении расчета на уровне HF/6-31G?

5.9. Работа №4. Расчет энергии активации реакции

Для расчета энергии активации реакции проводят ряд расчетов с фиксированным значением "координаты реакции". В качестве последней обычно выбирают расстояние между реакционными центрами реагента и субстрата. В результате определяется приблизительная геометрия и энергия переходного состояния, соответствующего максимуму на энергетической кривой. Точно переходное состояние рассчитывается путем задания параметра RUNTYP=SADPOINT и включения в исходный файл матрицы вторых производных в группу команд $HESS, которая располагается после группы $DATA. Для этого предварительно нужно рассчитать гессиан с приблизительной геометрией переходного состояния, после чего группу $HESS, содержащую матрицу вторых производных, копируют из файла с расширением ".dat", который создается после расчета в директории "DATA" (GAMESS-2003).

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

Это справедливо и для достаточно простых реакций, особенно если рассматриваются сольватационные эффекты. Например, в случае реакции S_N2 на поверхности потенциальной энергии присутствуют два минимума и один максимум (рис. 5.5).

Р и с. 5.5. Типичный вид профиля поверхности
потенциальной энергии S_N2 реакции: R – реагенты, Р – продукты,
RC – комплекс реагентов, PC – комплекс продуктов,
E^# – энергия активации реакции, H – энтальпия реакции,
E_с^# – величина центрального барьера

5.10. Пример 1. Расчет энергии активации S_N2 реакции бромметана с анионом хлора

1) построить график зависимости энергии от расстояния между атомом хлора и атомом углерода молекулы бромметана;

2) вычислить энергию активации реакции.

Расчет начинают с оптимизации геометрии комплекса реагентов (бромметан + хлорид-ион) и комплекса продуктов (хлорметан + бромид-ион). Это позволяет выбрать координату реакции, установить пределы и шаг её изменения, что определяет количество расчетов. Поскольку каждый расчет связан с оптимизацией геометрии, то от их количества зависит как общее время расчетов, так и часто сама возможность их проведения. В целях экономии времени большинство точек на кривой профиля поверхности потенциальной энергии рассчитывают в минимальном базисе. С наилучшей точностью вычисляют геометрии и энергии только экстремальных точек – минимумов (комплексы, интермедиаты) и максимумов (переходные состояния), что позволяет получать более точное значение энергии активации реакции. Мощность современных персональных компьютеров позволяет в качестве простейшего расчета использовать неэмпирический метод Хартри – Фока в базисе RHF/6-31G или RHF/6-31G(d). Для уточнения геометрии наиболее часто применяют метод МР2 в расширенном базисе, содержащем поляризационные и диффузные функции, а далее рассчитывают точную электронную энергию с учетом возможно большей части корреляционной энергии, используя методы теории возмущений более высокого порядка или метод связанных кластеров. Критерием является максимальное использование возможностей доступного исследователю компьютера.

Для изучаемой реакции исходный файл для построения графика зависимости энергии от координаты реакции и грубого расчета переходного комплекса имеет следующий вид:

$CONTRL SCFTYP=RHF RUNTYP=OPTIMIZE NZVAR=12 MAXIT=30

ICHARG=-1 MULT=1 COORD=ZMT $END

$SYSTEM TIMLIM=525600 MEMORY=1000000 $END

$BASIS GBASIS=N31 NGAUSS=6 $END

$SCF DIRSCF=.TRUE. $END

$STATPT OPTTOL=0.0001 NSTEP=20 IFREEZ(1)=10 $END

$DATA

C1

C

H 1 1.07349 2 106.8988

H 1 1.07356 3 111.9125 2 116.7619

H 1 1.07354 4 111.9156 3 -126.5581

Cl 1 3.500 2 179.000 3 55.0878

$END
Изменяя 10-ю координату, соответствующую расстоянию Cl-С, получим зависимость энергии от этого расстояния (табл. 5.5, расчеты №1 и №11 проведены без "замораживания" 10-й координаты, расчет №6 – в режиме оптимизации "седловой" точки).

Максимум (активированный комплекс, расчет №10 в табл. 5.5) точно рассчитывается в режиме RUNTYP = SADPOINT после вычисления гессиана (RUNTYP = HESSIAN) для максимума энергии при расстоянии C-Cl 2.4 Å:

$SYSTEM TIMLIM=525600 MEMORY=1000000 $END

$BASIS GBASIS=N31 NGAUSS=6 $END

$SCF DIRSCF=.TRUE. $END

$STATPT OPTTOL=0.0001 NSTEP=50 $END

$DATA

C1

C

H 1 1.0609759 2 89.1290284

H 1 1.0610023 3 119.9733469 2 88.4992711 0

H 1 1.0609646 4 119.9913347 3 177.0172667 0

CL 1 2.4000000 2 179.9998538 3 -124.8683605 0

$END

ENERGY IS -3068.8896422881 E(NUC) IS 169.1338948864

1 1-2.49609972E-02-7.66687296E-05 1.12593783E-04

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18 4-7.21736664E-06-7.68326794E-07 7.68683166E-03

$END

C

BR 1 2.5529423

H 1 1.0608874 3 120.0199635 2 89.4609332 0

H 1 1.0609282 4 119.9989036 3 178.9210329 0

CL 1 2.4104548 2 179.9993319 3 12.3100219 0

TOTAL ENERGY = -3068.8896302158
В результате ход реакции можно представить также в виде схем (рис. 5.6) или графика (рис. 5.7). График наглядно показывает положение комплекса исходных реагентов (минимум справа), комплекса продуктов реакции (минимум слева) и переходного состояния (максимум энергии).

Расстояние C-Cl 3.1 Å (комплекс реагентов)
Расстояние C-Cl 2.41 Å (переходный комплекс)
Расстояние C-Cl 1.9 Å (комплекс продуктов)

Р и с. 5.6. Структуры комплексов реакции CH₃Br + Cl^-

Энергия активации обратной реакции (хлористый метил + бромид-ион): 24.0 кДж/моль.

Экспериментальное значение энергии активации бромэтана с хлоридом лития в ацетоне: 18 ккал/моль. Почему так велико отличие расчета от эксперимента? Насколько велико различие энергии комплекса реагентов и изолированных реагентов и почему?

Р и с. 5.7. Зависимость энергии от координаты реакции CH₃Br + Cl^-

5.11. Пример 2. Расчет энергии активации S_N2 реакции бромметана с анионом хлора

1) построить график зависимости энергии от расстояния между ионом хлора и атомом углерода молекулы бромметана;

2) вычислить энергию активации реакции;

3) сравнить результаты расчета гидратированной и негидратированной систем.

Комплекс реагентов можно представить в виде следующей схемы:

Аналогично предыдущему примеру, сначала проводят расчет исходных реагентов, продуктов реакции и соответствующих комплексов (минимумов энергии) – оптимизируют их геометрии (RUNTYP = OPTIMIZE). Все расчеты обязательно нужно проводить в одном выбранном базисе и на одном и том же теоретическом уровне. Для возможности сравнения результатов с предыдущим примером будем проводить вычисления на уровне RHF/6-31G. Необходимо отметить, что уточнение энергетических параметров с использованием энергии нулевых колебаний и термохимических данных, получаемых при расчете гессиана (RUNTYP = HESSIAN), для уровня RHF/6-31G и исследуемой модели не приведет к лучшим результатам.

Изменяя координату, соответствующую расстоянию Cl-С, получают зависимость энергии от этого расстояния (табл. 5.6).

Далее рассчитывают геометрию и энергию переходного комплекса. Вычисленную приблизительно геометрию переходного состояния используют для точного поиска седловой точки, для чего предварительно рассчитывают гессиан (аналогично предыдущему примеру). Полученная зависимость энергии от расстояния C-Cl представлена в табл. 5.6 и на рис. 5.8. Ход реакции можно представить также в виде схем (рис. 5.9).

Величина центрального барьера (энергия "активации") составляет 22.4 кДж/моль, что меньше соответствующей величины для реакции негидратированных реагентов (25.9 кДж/моль). Из табл. 5.7 видно, что гидратация даже одной молекулой воды уменьшает разность энергий реагентов (продуктов) и их комплекса. Отсюда можно сделать вывод, что для расчета энергии активации даже в минимальном базисе необходимо учитывать сольватационные эффекты.
Таблица 5.6
Зависимость энергии от координаты реакции H2O^.CH₃Br + H2O^.Cl^-

№	Координата реакции (расстояние C-Сl), Å	Длина связи С-Br, Å	Энергия, а.е.	ΔE, кДж/моль
1 2 3 4 5 6 7 8 9 10 11 12 13 14	3.20 3.1010 min 3.00 2.80 2.60 2.50 2.44 2.4205 max 2.40 2.30 2.20 2.00 1.9463 min 1.90	2.0731 2.0843 2.0998 2.1495 2.2626 2.3952 2.50447 2.5419 2.5789 1.0627 2.8969 3.1282 3.1787 3.2255	-3220.9005838613 -3220,9007238100 -3220.9005395213 -3220.8987789453 -3220.8948362268 -3220.8927915851 -3220.8921985374 -3220.8921570042 -3220.8921998765 -3220.8935132216 -3220.8960894431 -3220.9011499917 -3220,9015454701 -3220.9011819946	2,52 2,16 2,64 7,26 17,61 22,98 24,54 24,64 24,53 21,08 14,32 1,04 0,00 0,95

Р и с. 5.8. Зависимость энергии от координаты реакции

H₂O^.CH₃Br + H₂O^.Cl^-
Таблица 5.7
Разность энергий реакций H2O^.CH₃Br + H2O^.Cl^- и CH₃Br + Cl^-

Частицы	Энергия, а.е.		Разность энергии комплекса и энергии реагентов/продуктов, а.е.
	Реагенты (продукты)	Комплекс реагентов (продуктов)
H2O.BrCH3+H2O.Cl-	-3220,885323699	-3220,900723810	-0,0154001104
BrCH3+Cl-	-3068,880630744	-3068,899488014	-0,0188572703
H2O.ClCH3+H2O.Br-	-3220,880411485	-3220,901545470	-0,0211339843
ClCH3+Br-	-3068,876507942	-3068,898768605	-0,0222606625

-535,5345421736 а.е.

-2685,3507815260 а.е.

Р и с. 5.9 (начало). Структура исходных реагентов и продуктов реакции
H2O^.CH₃Br + H2O^.Cl^-

-575,0482347388 а.е.

-2645,8321767470 а.е.