Gaussian: поиск TS и IRC на сложных PES (нужна методология)

[Metalian]

Считаю в гауссиане и потребовалось посчитать некоторую реакцию из металлокомплексного катализа. Судя по литературе, в подобных системах (да, реагент чуть другой и катализатор чуть другой, но можно для начала взять реагент, что был в статьях, и попробовать перенести на наш катализатор). Интересуют меня две стадии, собственно, A->B->C, и также есть вариант в одну стадию A->C. Экспериментально можно посмотреть только состав продуктов (есть ещё и побочная реакция, но с воспроизведением её механизма из литературы всё получше, и барьер там большой, так что можно пока забить), in situ ничего не посмотришь.

После долгих бесплодных мучений, а также последующего воспроизведения литературных данных я пришёл к следующему:
TS (A->B) и TS (A->C) лежат весьма близко по своей структуре и энергиям, и зачастую поиск TS (A->B) приводит к TS (A->C). Также иногда он приводит к TS (B->C) (и оно тоже не сильно далеко), что подтверждается IRC. Считать приходится, используя TS,CalcAll, зачастую ещё MaxStep уменьшать до 3-5 где-нибудь. QST2 вылетает почти всегда, QST3 тоже часто вылетает, возможно, в том числе и потому, что ротамеры и прочие конформеры свои чёрные дела делают (по меньшей мере сколько-то было случаев, когда спуск IRC приводит не к тому конформеру, с которого делался скан). Уменьшение MaxStep иногда помогает (видно, что первый шаг слишком большой), но чаще всего нет. Многочисленные сканы (между A, B и C), чаще всего приводят, в конечном итоге (после нахождения приемлемых приближений для поиска TS) к TS (A->C), реже к TS (B->C). Ну иногда может и не быть приемлемых вариантов вообще, такое тоже бывает, или скан уходит куда-нибудь не туда, ну да ладно. Также временами находится TS, похожее на правду (у меня несколько разных моделей катализаторов, поэтому разные случаи получаются), однако с мнимой частотой меньше 100 см^-1, и IRC сделать не получается (вылет на первом шаге или же останов после одного из шагов в каждую сторону). Ни настройки predictor,corrector, ни уменьшение шага с дефолтных 10 до 1-2, даже в сочетании с CalcAll не помогает.

Соответственно вопросы:
1. До каких значений целесообразно уменьшать MaxStep для поиска TS и какими ещё настройками можно поиграть, чтобы поиск TS из стартовой геометрии пошёл-таки в нужную сторону? При уменьшении с дефолтных 30 до 3-5 срабатывало иногда.
2. Есть ли какие-то способы успешно сделать IRC для случая, если единственная мнимая частота меньше 100 см^-1? С какими настройками поиграть?

[amge]

Я не знаю гауссиана, буду что-то говорить из общих соображений.

1. До каких значений целесообразно уменьшать MaxStep для поиска TS и какими ещё настройками можно поиграть, чтобы поиск TS из стартовой геометрии пошёл-таки в нужную сторону? При уменьшении с дефолтных 30 до 3-5 срабатывало иногда.

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

2. Есть ли какие-то способы успешно сделать IRC для случая, если единственная мнимая частота меньше 100 см^-1? С какими настройками поиграть?

Странно, что 100 см^-1. Вот меньше 10 - там да, проблемы.
Главная настройка - это, конечно, шаг IRC. У него размеренность - (корень квадратный из массы) * на длину, т.е. зависит от массы. Соответственно, если в ПС по мнимой моде перемещается легкий атом (например, водород), то шаг надо делать маленький, тяжелый - наоборот. Если IRC вылетело, сделав только 1-2 шага, то иногда помогает перезапустить его из последней геометрии с тем же гессианом. Ну и всегда можно с места обрыва IRC пустить обычную оптимизацию, хотя она негарантированно приведет в тот же минимум, что и IRC.

[Metalian]

Поиск ПС - черная магия, и гарантированно приводящих к успеху рецептов не существует.

Это понятно, и вроде бы опыт немалый, однако тут застрял.

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

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

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

Именно кривую поиска TS с некоторого начального приближения с подходящей мнимой частотой?

Кстати, заметил тут вот какую фигню: во всех расчётах, где удаётся из TS, которое вроде бы (A->B), спуститься не в C, а в B, кривая спуска имеет точки перегиба. Прикладываю одну из таких кривых.

paper-IRC-from-TS4-Def2SVP.log-graph.png

Если посмотреть на сами структуры, то видно, что сначала идёт движение по одной координате (да и по виду колебания с мнимой частотой в TS трудно предположить, что оно приведёт, куда надо), а потом вдруг движется по другой (и если посмотреть гессианы в точках, где идёт изменение направления, то там уже другая мнимая частота, которая какое-то время растёт по мере спуска, и она уже соответствует образованию интермедиата B).

2. Есть ли какие-то способы успешно сделать IRC для случая, если единственная мнимая частота меньше 100 см^-1? С какими настройками поиграть?

Странно, что 100 см^-1. Вот меньше 10 - там да, проблемы.

Тем не менее, Гауссиан что-то упорно не хочет делать спуск из такого TS.

Главная настройка - это, конечно, шаг IRC. У него размеренность - (корень квадратный из массы) * на длину, т.е. зависит от массы. Соответственно, если в ПС по мнимой моде перемещается легкий атом (например, водород), то шаг надо делать маленький, тяжелый - наоборот.

В тех стадиях, где перемещается водород, с дефолтными настройками всё хорошо, никаких проблем нет. Там, где речь идёт об образовании и/или разрыве связей C-C, тоже использовал дефолтный шаг. При слишком большом шаге велик шанс проскочить мимо вот этой вот хитрой области, где меняется направление движения, и сразу уйти в C.

Если IRC вылетело, сделав только 1-2 шага, то иногда помогает перезапустить его из последней геометрии с тем же гессианом.

Ну тут всякие опции с пересчётом гессиана через сколько-то шагов, вплоть до CalcAll. Но это если вылетело не на первом шаге, но при этом далеко от минимума.

Ну и всегда можно с места обрыва IRC пустить обычную оптимизацию, хотя она негарантированно приведет в тот же минимум, что и IRC.

Это точно, бывает, что оптимизация очень далеко уводит, и можно провраться.

[alex7851]

По какой конкретно причине вылетает IRC?
Из общих соображений - могу посоветовать попробовать еще Nudged Elastic Band (NEB), реализованный, например, в ORCA.

[Metalian]

По какой конкретно причине вылетает IRC?

Бывает так: Maximum number of corrector steps exceded. При этом опции стоят irc=(maxpoints=100,maxcycle=100,rcfc,recalcfC=(predictor=5,corrector=5)). И судя по всему, ни одного шага так и не сделано.

Бывает так: PES minimum detected on this side of the pathway - делает по шагу в каждую сторону и вылетает.

Существенное уменьшение шага не помогает.

Из общих соображений - могу посоветовать попробовать еще Nudged Elastic Band (NEB), реализованный, например, в ORCA.

Думал об этом. Может, тогда уж засунуть в периодику и посчитать NEB в привычном Quantum Espresso?

Только вот одна проблема с NEB: структуры TS он, в общем-то, обычно не даёт, её всё равно потом придётся дожимать. И не совсем понятно, как он обработает такую вот последовательно-параллельную реакцию (раньше доводилось NEBом некоторые элементарные реакции на поверхностях пытаться считать).

[amge]

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

Именно кривую поиска TS с некоторого начального приближения с подходящей мнимой частотой?

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

Кстати, заметил тут вот какую фигню: во всех расчётах, где удаётся из TS, которое вроде бы (A->B), спуститься не в C, а в B, кривая спуска имеет точки перегиба.

Это скрытые интермедиаты/переходные состояния. Иногда там бывают интересные/важные структуры. Обычно, если их обсуждают, то вместе с графиком энергии приводят график IRC-градиента, на котором в этих точках будут минимумы.

[amge]

Только вот одна проблема с NEB: структуры TS он, в общем-то, обычно не даёт, её всё равно потом придётся дожимать. И не совсем понятно, как он обработает такую вот последовательно-параллельную реакцию (раньше доводилось NEBом некоторые элементарные реакции на поверхностях пытаться считать).

У меня тоже сложилось впечатление, что NEB, QST2 и прочие подобные надежно работают только на простых реакциях, а на сложных ППЭ - фиг. Кроме того, для них нужно точно знать реактант и продукт, и не абы какие, а со строго соответствующей нумерацией атомов, что далеко не всегда просто. Гарантированный способ получить такие реактант и продукт - найти ПС и спуститься с него в обе стороны по IRC.

[alex7851]

Бывает так: Maximum number of corrector steps exceded. При этом опции стоят irc=(maxpoints=100,maxcycle=100,rcfc,recalcfC=(predictor=5,corrector=5)).

Попробуйте irc=(LQA, ...).

[amge]

По какой конкретно причине вылетает IRC?

Бывает так: Maximum number of corrector steps exceded. При этом опции стоят irc=(maxpoints=100,maxcycle=100,rcfc,recalcfC=(predictor=5,corrector=5)). И судя по всему, ни одного шага так и не сделано.

Бывает так: PES minimum detected on this side of the pathway - делает по шагу в каждую сторону и вылетает.

Существенное уменьшение шага не помогает.

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

[Metalian]

Это скрытые интермедиаты/переходные состояния. Иногда там бывают интересные/важные структуры. Обычно, если их обсуждают, то вместе с графиком энергии приводят график IRC-градиента, на котором в этих точках будут минимумы.

Спасибо! Прочитал. Но вот как это реализовать, не совсем понял (к сожалению, математики и физики в нужном объёме не знаю, чтобы понять, как это реализовать практически).

У меня тоже сложилось впечатление, что NEB, QST2 и прочие подобные надежно работают только на простых реакциях, а на сложных ППЭ - фиг. Кроме того, для них нужно точно знать реактант и продукт, и не абы какие, а со строго соответствующей нумерацией атомов, что далеко не всегда просто. Гарантированный способ получить такие реактант и продукт - найти ПС и спуститься с него в обе стороны по IRC.

Ну логично.
Строгая нумерация атомов - это понятно, и это как раз не очень сложно за ней уследить, и если надо, строчки поперемещать для единообразия.
Если бы в моём случае этот "гарантированный" способ сработал - не было бы этой темы.

Попробуйте irc=(LQA, ...).

Всё то же самое. По одному шагу и останов "по достижению минимума" на 1-м же шаге в каждую из сторон. Увеличить шаг?

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

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

[amge]

Это скрытые интермедиаты/переходные состояния. Иногда там бывают интересные/важные структуры. Обычно, если их обсуждают, то вместе с графиком энергии приводят график IRC-градиента, на котором в этих точках будут минимумы.

Спасибо! Прочитал. Но вот как это реализовать, не совсем понял (к сожалению, математики и физики в нужном объёме не знаю, чтобы понять, как это реализовать практически).

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

[Metalian]

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

Градиент на каждом шаге - это эта величина?
Maximum DWI gradient std dev = 0.022358046 at pt 64
А по абсциссе откладывать координату реакции, т.е. отсюда её брать?
NET REACTION COORDINATE UP TO THIS POINT = 0.24015

[amge]

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

Градиент в каждой точке IRC - это, по сути, наклон на графике энергии. Где он там в логе - не знаю, сам ни разу не считал гауссианом. Но точно знаю, что где-то есть, т.к. в одной из статей соавтор, считающий гауссианом, приводил такие графики. Судя по его блогу, там фигурируют слова "RMS Gradient Norm"

[Metalian]

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

Ну да, похоже на то (ну т.е. CHANGE IN THE REACTION COORDINATE чуть-чуть пляшет, но, думаю, это не настолько критично даже при большом количестве шагов, но да, наверное, проще использовать номер точки, хотя если скрипт сделать, то без разницы).

Градиент в каждой точке IRC - это, по сути, наклон на графике энергии. Где он там в логе - не знаю, сам ни разу не считал гауссианом. Но точно знаю, что где-то есть, т.к. в одной из статей соавтор, считающий гауссианом, приводил такие графики. Судя по его блогу, там фигурируют слова "RMS Gradient Norm"

Вижу такое в выхлопе:
Cartesian Forces: Max 0.003439657 RMS 0.000474061
Иного по шаблону "RMS" не находится.
Оно? Т.е. то, что после RMS?

[amge]

Вижу такое в выхлопе:
Cartesian Forces: Max 0.003439657 RMS 0.000474061
Иного по шаблону "RMS" не находится.
Оно? Т.е. то, что после RMS?

Похоже. Можно проверить, если минимумы будут там, где у энергии перегибы, то, скорее всего, оно.

PS Неожиданно обнаружил у себя свой же скрипт fchk2xyz, который вычисляет градиенты из содержимого поля 'IRC point 1 Gradient at each geome'. Насколько правильно - не знаю. Наверное, проверял, но уже не помню.

[Metalian]

Ну если я всё правильно извлёк, то получается примерно так.
По оси абсцисс - координата реакции (не номер точки), по оси ординат - Cartesian Forces: RMS. Сглаживать не стал.

Это для реакции A->B:

A-B.png

Это для реакции A->C:

A-C.png

Но это для разных ротамеров предреакционного комплекса, и с поиском переходного состояния B -> C возникли проблемы. Главное здесь, что A-B удалось.

Однако когда я взял найденное "удачное" TS A->B, крутанул на 180 градусов вокруг связи металл - реагент (речь о металлокомплексном катализе), после оптимизации TS почти не изменилось (мнимая частота и геометрия очень близки к тому, что было), спуск привёл-таки к A и C, и без всяких скрытых минимумов.

A-C-rot180.png

Есть идеи? Брать точки, соответствующие минимумам градиента и искать там TS?

[Metalian]

Попробовал для реакции A-B (положительная координата реакции ведёт к A, отрицательная к B) оптимизацию из точки, соответствующей промежуточному пику - получил B.
Также попробовал найти TS из промежуточного минимума градиента - получилось TS c мнимой частотой 110 см_-1, что-то в ней явно есть, но попытка сделать IRC ожидаемо ни к чему не провела: в обе стороны на первой же точке программа говорит о достижении минимума.
Похожее TS недавно уже было у меня, и не помогло ничего: ни увеличение или уменьшение шага, ни CalcAll, ни Cartesian, ни Tight, ни Phase, ни LQA.
Есть идеи, как ещё можно что-то с ним сделать?

[amge]

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

[Metalian]

Ладно, а что делать с запуском IRC из TS, в котором мнимая частота в районе 100 см^-1, чтобы это IRC всё же пошло, а не остановилось на первом шаге?

Другой прикол: из TS (получено с CalcAll, так что вряд ли есть смысл его "дожимать") в обе стороны получился один и тот же интермедиат (а должен был в одну из сторон предреакционный комплекс). Пересчитать IRC c CalcAll?

Ну и третий: путь из TS к предреакционному комплексу пошёл в его сторону, но не дошёл, и свернул к интермедиату ещё одного варианта реакции (помимо A->C и A->B->C есть ещё A->D, но до этого пути этих реакций не пересекались). Ну это ладно, переживём, посчитаем отдельно...

[amge]

Мог бы посоветовать перейти на Природу - у нее с IRC как правило проблем нет, но не буду, по двум причинам.
Во-первых про проблемы с IRC у гауссиана тоже раньше не слышал.
Во-вторых, здесь, видимо, сложная ППЭ с близко расположенными ПС. Сменишь хоть немного уровень теории, и все изменится до неузнаваемости.

Так что я больше ничего не могу посоветовать.