ПРИРОДА. Знакомство с основами расчёта на примере конформационного анализа NH2-OH

На примере конформационного анализа гидроксиламина (NH₂-OH) показывается последовательность квантовохимического расчёта. Приведены входные данные для квантовохимического пакета ПРИРОДА, версия 06.

1. Оптимизация структуры соединений

На первом этапе производится поиск структуры, отвечающей локальному минимуму на поверхности потенциальной энергии (ППЭ) вблизи начально заданной структуры (пробной структуры). В итоге будут получены геометрические параметры одного из конформеров NH₂-OH и значение полной энергии в атомных единицах (хартри).

$system memory=512 disk=10 path=. $end
$control theory=dft task=optimize basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.026
1 2 1.026 1 105.116
8 2 1.464 3 102.846 1 107.322
1 4 0.972 2 101.187 1 -125.643
$end

2. Решение колебательной задачи

Для проверки соответствия оптимизированной структуры точке локального минимума на ППЭ производится расчёт гессиана (матрицы вторых производных энергии по координатам ядер). Если в диагонализированной матрице отсутствуют мнимые силовые постоянные (частоты колебаний), значит структура отвечает локальному минимуму на ППЭ. Фактически также происходит расчёт ИК-спектра. Ввиду малой затратности методами статистической термодинамики (ТД) также производится расчёт энтропии (S) и термохимических поправок (относительно полной энергии E) к энтальпии (H), внутренней энергии (U, за вычетом энергии ядер) и потенциала Гиббса (G) в ккал/моль. Имеет смысл сразу расчитывать поправки для всех температур, которые гипотетически потребуются в будущем (т.к. повторный расчёт занимает несоизмеримо больше времени, чем расчёт ТД поправок).

$system memory=512 disk=10 path=. $end
$control theory=dft task=hessian basis=»3z.bas» $end
$dft functional=PBE $end
$thermo temperature=0.0, 10.0, 20.0, 298.15 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.026
1 2 1.026 1 105.092
8 2 1.464 3 102.886 1 107.360
1 4 0.972 2 101.187 1 -125.604
$end

3. Релаксированное сканирование геометрических параметров

Для поиска оставшихся конформеров необходимо произвести сканирование угла поворота группы -OH относительно -NH₂. Значение 3 в параметре fix говорит о том, что изменяемый геометрический параметр является двугранным (торсионным) углом. В конце выходного файла будет приведён график сканирования в координатах x — E.

$system memory=512 disk=10 path=. $end
$control theory=dft task=scan basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 fix=3,5,4,2,1 value=-125.6, 54.4 points=10 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.026
1 2 1.026 1 105.116
8 2 1.464 3 102.846 1 107.322
1 4 0.972 2 101.187 1 -125.643
$end

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

$system memory=512 disk=10 path=. $end
$control theory=dft task=optimize basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.027
1 2 1.027 1 106.838
8 2 1.449 3 106.700 1 114.122
1 4 0.978 2 107.367 1 54.400
$end

По аналогии с расчётом первого локального минимума (который является и глобальным минимумом на ППЭ вращения группы -OH), нужно решить колебательную задачу для второго минимума, проверить, что оптимизированная структура действительно отвечает минимуму на ППЭ.

$system memory=512 disk=10 path=. $end
$control theory=dft task=hessian basis=»3z.bas» $end
$dft functional=PBE $end
$thermo temperature=0.0, 10.0, 20.0, 298.15 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.027
1 2 1.027 1 107.071
8 2 1.448 3 106.900 1 114.275
1 4 0.978 2 107.367 1 57.295
$end

4. Оптимизация переходного состояния

Для расчёта переходного состояния реакции превращения одного конформера в другой нужно извлечь прибжённую структуру, отвечающую максимуму полной энергии из файла релаксированного сканирования, и совершить с ней ряд действий. Первое — нужно произвести расчёт пробного гессиана для взятой структуры. В выходном файле необходимые данные будут заключатся в секции $energy … $end. В этой секции из каждой строки необходимо удалить префикс ‘MOL>’

$system memory=512 disk=10 path=. $end
$control theory=dft task=hessian basis=»3z.bas» $end
$dft functional=PBE $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.029
1 2 1.029 1 103.549
8 2 1.476 3 102.089 1 110.561
1 4 0.973 2 105.825 1 -5.600
$end

Второе, это непосредственно поиск переходного состояния (ПС). В качестве входных данных используется структура, отвечающая максимуму на ППЭ и рассчитанный для неё гессиан.

$system memory=512 disk=10 path=. $end
$control theory=dft task=optimize basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 saddle=1 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.029
1 2 1.029 1 103.549
8 2 1.476 3 102.089 1 110.561
1 4 0.973 2 105.825 1 -5.600
$end
$energy
e=-1.316211396765e+02
d=-9.65246921e-02 -7.46838060e-01 8.79204087e-01
a=-1.92680688e+01
…
$end

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

$system memory=512 disk=10 path=. $end
$control theory=dft task=hessian basis=»3z.bas» $end
$dft functional=PBE $end
$thermo temperature=0.0, 10.0, 20.0, 298.15 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.029
1 2 1.029 1 103.490
8 2 1.478 3 101.706 1 110.176
1 4 0.973 2 105.493 1 -11.599
$end

5. Сканирование по внутренней координате реакции

Для окончательного убеждения в том, что найденное ПС лежит на ППЭ реакции превращения одного конформера в другой, т.е. связывает найденные локальные минимумы, нужно произвести сканирование вдоль внутренней координаты реакции (IRC). В качестве входных данных используется оптимизированная структура ПС и расчитанный для неё гессиан. Также во входном файле указывается направление спуска из ПС (back=0 или back=1 в группе $optimize), спуск производится по направлению наибольшей кривизны ППЭ.

$system memory=512 disk=10 path=. $end
$control theory=dft task=irc basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 back=0 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.029
1 2 1.029 1 103.490
8 2 1.478 3 101.706 1 110.176
1 4 0.973 2 105.493 1 -11.599
$end
$energy
e=-1.316210637441e+02
d=-9.47079233e-02 -7.53916266e-01 8.30176825e-01
a=-1.92994298e+01
…
$end

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

$system memory=512 disk=10 path=. $end
$control theory=dft task=irc basis=»3z.bas» $end
$dft functional=PBE $end
$optimize steps=100 back=1 $end

Можно рассчитать абсолютные и относительные (стандарт Si(CH₃)₄, абсолютные хим.сдвиги для стандарта указываются в группе $nmr) химические сдвиги ЯМР. Необходимо помнить, что относительные величины воспроизводятся квантовой химией значительно лучше, чем абсолютные. Поэтому лучшим параметром для сравнения с экспериментом будет положение сигналов относительно друг друга. Единичный расчёт хим.сдвигов не учитывает конформационную подвижность различных групп (усреднение и уширение сигналов ЯМР).

$system memory=512 disk=10 path=. $end
$control theory=dft task=nmr basis=»3z.bas» $end
$dft functional=PBE $end
$nmr standard(1)=31.3867 standard(6)=181.9781 $end

$molecule charge=0 mult=1 z-matrix
1
7 1 1.026
1 2 1.026 1 105.092
8 2 1.464 3 102.886 1 107.360
1 4 0.972 2 101.187 1 -125.604
$end

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

Вы можете скачать все расчёты одним файлом (171 KiB).

Тема на нашем форуме, посвящённая программе ПРИРОДА.

Панкратьев Евгений, 2008