Координата реакции

[eFo]

Уважаемые коллеги, может ли координата реакции иметь такой вид? Или всё-таки в области плато между IM и TS1-b нужно попытаться найти ещё одно TS (но там всё-таки идёт постоянное увеличение энергии)? Такая координата получилась анализом IRC...

[sanya1024]

Посмотрите, какая структура отвечает этому плечу. Каким внутренним координатам отвечает координата реакции? (очевидно, на разных участках главный вклад в коорд. реакции дают разные внутренние координаты).

[amge]

Красиво. По моему, все нормально. Между IM и TS1-b на плато может быть ПС, но может и не быть. Поискать, конечно, стоит. А также на плече сплошной линии между Pr1 и TS2. Если ПС не обнаружатся, но захочется пообсуждать "плечевые" структуры, то советую привести еще один график, на котором на ординате не энергия, а градиент IRC, примеры можно найти в блоге Henry Rzepa.

Хочу еще заметить, что шкала абсцисс немного двусмысленна. Например, на пунктирной линии различие между продуктами 90 - это, вероятно, количество шагов IRC умноженное на величину шага. Но его можно трактовать и как действительное различие между продуктами (в тех же единицах, корень момента инерции между ними, эту величину тоже можно вычислить, она, скорее всего, будет меньше, как длина ломаной линии vs. расстояние между ее концами). И даже два Pr2 с различными абсциссами не спасают: в результате двух реакций, приводящих к одинаковому продукту, нумерация атомов в нем может стать различной, и моменты инерции тоже. По моему, лучше на оси асцисс давать количество шагов IRC и отдельно указать величину шага.

[o-oxhem]

Правильная кривая IRC должна иметь форму колоколообразной кривой. Отклоненения являются или признаком недостаточной точности расчета, или существованием точки бифуркации ПМЭР , VRI-точки. Нужно исследовать эти фрагменты ППЭ подробнее.
Далее по рис. - такое изображение не совсем верно, поскольку внутренние координаты могут быть разными и откладывать их на одной оси не следует, лучше и правильнее, чтобы не было дублирования одного интермедиата на рис., а ось оставить без подписи и без размерностей.

[amge]

Меня давно интересовало, можно ли и если да, то как найти на кривой IRC точки бифуркации (VRI-точки), или хотя бы предположить их местоположение, но без привлечения других методов. И тут под впечатлением поста o-oxhem решил провести небольшое расследование. Оказалось, не я первый, и вывод такой: IRC сам по себе непригоден для поиска VRI. Тем не менее, захотелось убедится. Взял простейшую классическую систему с бифуркацией, перегруппуровку метоксильного радикала в гидроксиметиленовый, сделал IRC и посчитал гессианы в каждой точке. Результат неутешительный: VRI-точка (на синем рисунке она помечена кружочком) оказалась совсем не там, где я мог бы предположить. Такие упоминающиеся в литературе признаки VRI, как перегиб на кривой IRC и нарушение симметрии нифига не работают: до перегиба еще довольно далеко, а симметрия не нарушается, т.е. IRC приходит к точке бифуркации и тупо шпарит по гребню, пока не достигает второго переходного состояния (оно на пологом плече).

VRI.png

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

[sanya1024]

Я правильно понимаю, что белый рисунок -- что должно быть, а синий -- что получилось? Логично: исходно симметричная структура 1 через симметричное же переходное состояние 2... что делает? то ли доходит до опять-таки симметричной структуры 3, к-рая снова оказалась переходным состоянием, только по другой внутренней координате, то ли сразу за точкой бифуркации начинает сползать в одну из несимметричных ям 4 или 4'. А у Вас задача IRC остановилась в точке 3? Т.е., если теперь ее рестартнуть, выбрав новую моду с мнимой частотой, то IRC поползет уже к точкам 4 или 4', так? А если потом склеить обе кривые, откладывая по абсциссе не значение внутренней координаты, а номер шага IRC, то и получим кривую с плечом, как в исходном топике.

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

А вот почему из точки бифуркации IRC процедура не приводит сразу к точкам 4 или 4', а ведет по гребню? видимо, дело в алгоритме. Алгоритм говорит: следим за модой такой-то (той, что имела мнимую частоту в точке 2). А то, что у нас по ходу дела появилась новая мода с мнимой частотой -- это алгоритм не волнует. Ему за этой модой следить не велено Может, где-то есть более "умные" алгоритмы?

[amge]

А у Вас задача IRC остановилась в точке 3? Т.е., если теперь ее рестартнуть, выбрав новую моду с мнимой частотой, то IRC поползет уже к точкам 4 или 4', так?

У меня IRC прошла через точку 3 (это ПС, оно на пологом плече), но не остановилась на ней и дошла до конца (до 4 или 4'). Насколько я понял, остановится IRC на переходном состоянии или нет, зависит от случая. Так, то что приведено на синем рисунке - IRC с шагом 0.01 A.amu^1/2 - проскочило. С шагом 0.05 - остановилось. С шагом 0.2 - опять проскочило. Видимо, если IRC на очередном шаге случайно попадает почти точно на ПС, где все градиенты нулевые, то там и заканчивается, иначе - проскакивает.

А вот почему из точки бифуркации IRC процедура не приводит сразу к точкам 4 или 4', а ведет по гребню? видимо, дело в алгоритме. Алгоритм говорит: следим за модой такой-то (той, что имела мнимую частоту в точке 2). А то, что у нас по ходу дела появилась новая мода с мнимой частотой -- это алгоритм не волнует. Ему за этой модой следить не велено Может, где-то есть более "умные" алгоритмы?

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

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

[sanya1024]

Как-то давно столкнулась с забавной ситуацией. Делала relaxed scan по внутренней координате и IRC по координате реакции, к-рая была очень похожа на эту внутреннюю координату, но все-таки не совпадала с ней. А энергию на графике откладывала в зависимости от этой внутренней координаты (т.е., для IRC это была как бы проекция на эту внутреннюю координату). Вроде бы IRC должен выбирать точки с минимальной энергией при заданных ограничениях, но на графике получилось, что иногда кривая IRC проходит чуть выше, чем кривая relaxed scan. Я пыталась разобраться, в чем дело, но дальше утверждения "ну там ограничения разные, поэтому и кривые по-разному прошли" не продвинулась. Наоборот, мне кажется, что на пути IRC меньше ограничений, там все подчинено только минимизации энергии. А relaxed scan -- это минимизация энергии с фиксированной одной координатой. Лишнее ограничение -- меньше вариационной свободы -- выше энергия. В основном (кроме небольших участков) так оно и было. Но почему на этих небольших участках было наоборот -- загадка. Для численных погрешностей эффект великоват (где-то сотые доли ккал/моль).

[asp]

Вопрос в тему. Пытаюсь посчитать координату реакции ацилирования амина хлорангидридом карбоновой кислоты. Считаю в Природе.
1. Сначала оптимизирую геометрии амина и хлорангидрида отдельно.
2. Помещаю молекулы рядом (около 3 А) и оптимизирую геометрию. В результате молекулы сближаются, но существенного взаимодействия не происходит, получаем предреакционный комплекс (прошу поправить если неправильно называю).
3. Гессианы во всех случаях не содержат мнимых частот.
4. Вручную сближаю молекулы и при этом искажаю углы так, чтобы стало похоже на тетраэдрическое переходное состояние (взял в Марч Дж. "Органическая химия: реакции, механизмы и структура" т.2 М.:Мир 1987) и считаю гессиан. Получаю мнимую частоту.
5. Запускаю расчет седловой точки. Вот здесь получается непонятное мнимая частота исчезает. Так может быть? И значит ли это, что нарисованная структура не является переходным состоянием?
6. Поступаю по другому, беру структуру после нескольких шагов оптимизации при поиске седловой точки (пока мнимая частота не исчезла) и запускаю расчет IRC. Получаю в одну сторону расчет идет, а в другую нет. См. файл.
Хотелось бы услышать советы более опытных.

[EvgeniX]

В общем случае ничто не запрещает ППЭ иметь хитрую форму со своеобразными "тупичками". В таких случаях как раз можно оптимизировать ложное ПС, но IRC в обоих случаях будет идти в одну сторону.

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

[eFo]

5. Запускаю расчет седловой точки. Вот здесь получается непонятное мнимая частота исчезает. Так может быть? И значит ли это, что нарисованная структура не является переходным состоянием?

Может конечно, но это совсем не означает, что нарисованная структура не является ПС. У меня такое бывает, не туда оптимизация уходит по ППЭ. Тогда ставлю в $optimize trust=0.001, оптимизирую, смотрю в визуализаторе, пока не исчезнет отрицательная кривизна, беру ближайшую структуру, считаю гессиан и заново.. Обычно всё помогает.

[eFo]

Да, совсем замотался, всем откликнувшимся спасибо за советы и подсказки, много нового для себя узнал. Более подробное исследование этих участков не дало никаких результатов (в плане ПС или ИМ)..

[amge]

4. Вручную сближаю молекулы и при этом искажаю углы так, чтобы стало похоже на тетраэдрическое переходное состояние (взял в Марч Дж. "Органическая химия: реакции, механизмы и структура" т.2 М.:Мир 1987) и считаю гессиан. Получаю мнимую частоту.
5. Запускаю расчет седловой точки. Вот здесь получается непонятное мнимая частота исчезает. Так может быть? И значит ли это, что нарисованная структура не является переходным состоянием?

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

Попробовал посчитать Вашу систему. Оказывается, Вы уже практически нашли ПС. У меня оптимизация (с 1-й точкой Вашего IRC в качестве начального приближения) без проблем за несколько шагов дала истинное ПС. С него опять же без проблем проходит IRC в обе стороны. Ниже - инпут с опциями, которые я использовал, и оптимизированной геометрией ПС.

Код: Выделить всё

$system memory=512 disk=-2000 $end
$control
 task=hessian
 theory=DFT
 basis=L1                                             
 print=+molden
$end
$dft functional=PBE $end
$grid accuracy=1e-8 $end
$optimize
 tolerance=1e-5
 trust=0.5
 steps=300
 radius=0.4
 points=100
$end
$molecule 
 cartesian
   6    0.60398117   -0.30036094    1.59220471
   8    1.70406897   -0.55186099    1.98711511
  17   -1.12844443   -0.27399274    3.32680357
   6   -0.22944022   -1.08483796    0.61535183
   1   -0.27280112   -2.11815884    0.98297841
   1    0.27481157   -1.08316803   -0.36803533
   1   -1.25262621   -0.70148664    0.53369482
   7    0.21902091    1.34773386    1.51887827
   6    1.09141231    2.06729621    2.47376553
   6    0.23497122    1.91414222    0.15244655
   1   -0.73277693    1.24517626    1.96650498
   1    0.85278703    1.67993354    3.47400295
   1    2.14120225    1.85123334    2.23259027
   1    0.88746200    3.14624328    2.41493281
   1   -0.49199978    1.38442030   -0.47686226
   1    1.24319760    1.80399827   -0.27544776
   1   -0.02814075    2.98253635    0.18780420
$end

Если понадобится, выложу и остальные файлы.

[asp]

Спасибо всем за советы. Посчитал с настройками от amge получилось. Но все равно загадкой осталось почему была выбрана именно та структура. Методом научного тыка? Хочу для себя понять каков алгоритм поиска ПС. Система была взята для тренировки.

[eFo]

Спасибо всем за советы. Посчитал с настройками от amge получилось. Но все равно загадкой осталось почему была выбрана именно та структура. Методом научного тыка? Хочу для себя понять каков алгоритм поиска ПС. Система была взята для тренировки.

А что плохого, если вы с помощью химической логики предложили ПС? Ну, а в общем случае relaxed scan по ППЭ.

[asp]

Плохого ничего в методе научного тыка нет. Сам им пользуюсь. Но хочется более быстрого метода. Не могли бы привести примеры relaxed scan по ППЭ для Природы в декартовых координатах?

[amge]

Не могли бы привести примеры relaxed scan по ППЭ для Природы в декартовых координатах?

Ниже пример инпута скана для Вашей тестовой системы: сближение молекул хлорангидрида и амина. Замечу, что в Природе скан устроен очень удобно: начальное приближение можно задавать хоть декартовыми координатами, хоть z-матрицей, и нет нужды заботиться о том, чтобы сканируемый параметр присутствовал в явном виде в z-матрице или в списке делокализованных координат.

Код: Выделить всё

$system mem=512 disk=-2000 $end                                    
$control
 task=scan
 basis=L1
$end
$grid accur=1e-8 $end
$optimize
 fix=1,1,8,0,0 value=2.7,1.5 points=13
 steps=300 tol=1e-5 trust=0.5
$end
$molecule 
 cartesian
   6    0.34489490   -0.84044039    2.01682932
   8    1.48494258   -0.83342010    2.33999510
  17   -0.99438773   -0.54764264    3.30772586
   6   -0.27803099   -1.22433283    0.71102455
   1   -0.65149670   -2.25768708    0.80390895
   1    0.48488697   -1.17910742   -0.07821650
   1   -1.12410523   -0.56859964    0.48003967
   7    0.19155072    1.90085177    1.26962101
   6    1.19972545    2.41506963    2.19379455
   6    0.30041537    2.48880652   -0.06208078
   1   -0.74045660    2.08268552    1.65471779
   1    1.02312134    1.99702787    3.19646619
   1    2.19876352    2.08062111    1.86605863
   1    1.21963104    3.52468726    2.26775137
   1   -0.52040028    2.12649849   -0.70201472
   1    1.25113114    2.17158913   -0.52499995
   1    0.28256184    3.60161973   -0.07404263
$end

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

file480x480.gif

Кстати, у этой реакции есть продолжение: образование более стабильной структуры, в которой HCl координирован не по азоту, а по кислороду. Интуитивно кажется, что это должен быть легкий процесс, но найти его ПС у меня с ходу не получилось...

[asp]

Провел сканирование, вопрос по параметрам fix=1,1,8,0,0 value=2.7,1.5
частично разобрался 1, 8 это первый и восьмой атомы, 2.7 начальная дистанция между ними, а что значат 1 перед 1, 8 и 1.5 не понял.
И еще в какой программе Вы делали анимацию?

[amge]

fix=1,1,8,0,0 Первая 1 означает, что сканируемый параметр - длина "связи". На ее месте 2 - валентный угол, 3 - диэдральный.
value=2.7,1.5 Начальное и конечное значения сканируемого параметра (углы - в градусах)
Как я делаю анимации

[asp]

Благодарю, за подробный ответ и терпение с нами, неучами.