Воспроизведение расчетов, описанных в статье

[Shadin]

Здравствуйте!
Пытаюсь повторить расчеты, которые были выполнены в статье: J.Phys.Chem.A 2009 p.2578 (Theoretical Study on the Mechanism and Diastereoselectivity of NaBH4 Reduction; Yasumitsu Suzuki, Daisuke Kaneno and Shuji Tomoda)
Считаю в FireFly структуру, приведенную на рис 6. Геометрию переходного состояния брала из Support Information.
Метод DFT (B3LYP), базис 6-31G(d)+.
С ПС все нормально, делаю IRC в обе стороны, тоже все нормально, затем оптимизирую продукты и реагенты. С реагентами проблема: во-первых BH4-группа уезжает, во-вторых мнимые частоты. Пробовала добавлять nonvdw и уменьшать opttol, ставить INTTYP=HONDO, ICUT=11, ITOL=30 в $CONTRL и FDIFF=.FALSE. в $SCF, но ничего не меняется.

Вопрос, что не так?

Спасибо огромное!
Извините, если ответ на поверхности!

[VTur]

Попробуйте для реагентов стартовать с минимума с помощью CONOPT в $statpt. Потом сравните переходные состояния.
Так непонятно, но может быть есть конкурирующий путь. Можно посоветовать еще уменьшить OPTTOL=0.0001, ну и другие стандартные вещи, типа учитывать больше интегралов.

[sanya1024]

При численном расчете гессиана могут получаться мнимые частоты на ровном месте. Поставьте $FORCE NVIB=2 PURIFY=.t. $END -- может помочь. Здесь на форуме уже объяснялось и сама процедура численного расчета гессиана, и значение этих параметров. Если у Вас сохранились файлы IRCDATA или *.irc, скопируйте из них (не из *.out файлов!) группу $VIB целиком в конец инпута, убедитесь, что в конце стоит $END (если нет, то поставьте). Тогда Вам не придется пересчитывать весь гессиан заново, а только досчитать недостающие 3N смещений.
Если все-таки не все мнимые частоты исчезли, слегка крутаните соответствующие группы (в Ваших примерах все мнимые моды -- внутренние вращения) руками в ChemCraft-е и переоптимизируйте.

Теперь "уезжающая" BH₄^- группа. Предполагается, что BH₄^- должен сесть водородом на карбонильный углерод -- как на рис. 7 в статье. Я вижу, что путь IRC приводит к такой структуре, а дальнейшая ее оптимизация приводит к "уезжанию". Значит, так тому и бывать: "уехавшая" структура имеет более низкую энергию, чем та, что в статье (см. выдачи). А кто сказал, что потенциальная поверхность этой системы должна быть как в учебнике: две долины и седло? там и долин может быть побольше, и долина может быть не долиной, а каким-то подобием плато, с к-рого есть путь в долину пониже.

[Shadin]

Попробуйте для реагентов стартовать с минимума с помощью CONOPT в $statpt. Потом сравните переходные состояния.
Так непонятно, но может быть есть конкурирующий путь. Можно посоветовать еще уменьшить OPTTOL=0.0001, ну и другие стандартные вещи, типа учитывать больше интегралов.

При численном расчете гессиана могут получаться мнимые частоты на ровном месте. Поставьте $FORCE NVIB=2 PURIFY=.t. $END -- может помочь. Здесь на форуме уже объяснялось и сама процедура численного расчета гессиана, и значение этих параметров. Если у Вас сохранились файлы IRCDATA или *.irc, скопируйте из них (не из *.out файлов!) группу $VIB целиком в конец инпута, убедитесь, что в конце стоит $END (если нет, то поставьте). Тогда Вам не придется пересчитывать весь гессиан заново, а только досчитать недостающие 3N смещений.
Если все-таки не все мнимые частоты исчезли, слегка крутаните соответствующие группы (в Ваших примерах все мнимые моды -- внутренние вращения) руками в ChemCraft-е и переоптимизируйте.

Спасибо большое за помощь, все это я попробую, сделаю! А про NVIB и PURIFY - совсем забыла про эти опции, простите.

Теперь "уезжающая" BH₄^- группа. Предполагается, что BH₄^- должен сесть водородом на карбонильный углерод -- как на рис. 7 в статье. Я вижу, что путь IRC приводит к такой структуре, а дальнейшая ее оптимизация приводит к "уезжанию". Значит, так тому и бывать: "уехавшая" структура имеет более низкую энергию, чем та, что в статье (см. выдачи). А кто сказал, что потенциальная поверхность этой системы должна быть как в учебнике: две долины и седло? там и долин может быть побольше, и долина может быть не долиной, а каким-то подобием плато, с к-рого есть путь в долину пониже.

Все это я понимаю. Однако авторы статьи пишут All geometries
of reactants, products, and transition state structures for the
reaction of NaBH4 with acetone and several substituted cyclohexanones
were optimized by using DFT calculations with
B3LYP hybrid functional and 6-31+G(d) basis set ... Stationary points were fully optimized and characterized
by vibrational frequency calculations
И так как структуру на рисунке 7 они называют реагентом, то видимо это оптимизированная геометрия? нет?

Также не очень понятно, почему авторы оптимизируют геометрии DFT (B3LYP), а энергию структур считают МР2, и почему не использую базис с диффузными функциями для водорода? Или эти вопросы нужно задавать не здесь, а авторам? Просто хочу полностью разобраться в приведенных в статье расчетах.

[VTur]

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

Что бы я сделал
- один совет я Вам дал
- резко увеличил бы базис
- если есть возможность, перешел бы на ГАУССИАН, там есть процедура QST2, которая стартует от известных реагентов и продуктов, сама располагает переходное состояние и строит путь реакции в заданном числе точек

[sanya1024]

Все это я понимаю. Однако авторы статьи пишут All geometries
of reactants, products, and transition state structures for the
reaction of NaBH4 with acetone and several substituted cyclohexanones
were optimized by using DFT calculations with
B3LYP hybrid functional and 6-31+G(d) basis set ... Stationary points were fully optimized and characterized
by vibrational frequency calculations
И так как структуру на рисунке 7 они называют реагентом, то видимо это оптимизированная геометрия? нет?

Также не очень понятно, почему авторы оптимизируют геометрии DFT (B3LYP), а энергию структур считают МР2, и почему не использую базис с диффузными функциями для водорода? Или эти вопросы нужно задавать не здесь, а авторам? Просто хочу полностью разобраться в приведенных в статье расчетах.

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

Авторы использовали Гауссиан, да еще версии 03, а это совершенно не способствует активным поискам на ППЭ: посчиталось что-то -- и ладно. Скорее всего, в долине реагентов картинка следующая:

reagent_valley.jpg

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

Почему авторы все оптимизации делали в B3LYP, а энергии считают в MP2 -- понятно. Оптимизация в MP2 -- дело дорогое и долгое даже в FireFly, что уж говорить про не приспособленный к этому Гауссиан, тем более версии 03. По опыту моему и моих коллег, расчеты DFT и MP2 провираются в противоположные стороны при оценке межмолекулярных взаимодействий, так что истина на самом деле где-то посередине. Я не исключаю, что методом CCSD(T) ("золотой стандарт" для межмол. взаимодействий -- в т.ч. и по затратам) или эмпирическим, но калиброванным по CCSD(T) методом DFT-D (дисперсионная поправка Гримме) получится, что все-таки минимум, полученный авторами статьи, чуточку глубже. Но это надо проверять. В Природе есть параллельный CCSD (без T), в GAMESS-US есть параллельный CCSD(T) для основного состояния. В Гауссиане, в принципе, он тоже есть, но боюсь, что такой же медленный и затратный. Короче, если перепроверять, то заведомо более высокоуровневым методом CCSD(T).

Насчет базисов тоже понятно: в Гауссиане не очень-то развернешься с большими базисами: очень быстро задача станет совершенно неподъемной.

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

[amge]

Пытаюсь повторить расчеты

Если Вы хотите воспроизвести чужую работу, постарайтесь ничего не менять по сравнению со статьей: ни базис, ни метод. А программу сменить можно

Относительно последнего тоже надо проявлять осторожность, если речь идет о восвроизведении результатов, полученных гауссианом. Известно, например, что гауссиан использует не вполне стандартную версию B3LYP (в последних версиях GAMESS US это B3LYP1). Есть ли аналогичная модификация в FireFly - не знаю.
Потом, результаты, полученные гауссианом даже одной и той же версии, но разных ревизий, могут заметно отличаться.

Так что, если воспроизводимость в общих чертах есть, то, наверное, можно этим и ограничиться.

[sanya1024]

Из мануала FireFly:

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

hybrid functionals
= B3LYP1 B3LYP as implemented in NWCHEM and GAUSSIAN 98, using VWN formula 1 RPA correlation
= B3LYP5 B3LYP as implemented in GAMESS (US), using VWN formula 5 correlation
= B3LYP either B3LYP1 or B3LYP5 depending on the value of B3LYP keyword in the $DFT group, see later.

Ну, и в GAMESS аналогично.
Про смену программы (точнее, оптимизатора) я сказала не случайно: если результат численно стабильный, а минимум -- достаточно глубокий, то смена оптимизатора не повлияет. Очевидно, в данном случае как раз авторы статьи попали в неглубокий минимум, а оптимизатор FF нашел минимум поглубже.

[Shadin]

Здравствуйте! Не писала, так как занималась расчетами, сейчас попытаюсь представить результаты, ну и, конечно же, у меня опять вопросы

Если Вы хотите воспроизвести чужую работу, постарайтесь ничего не менять. А программу сменить можно

К GAUSSIAN ХХ у меня вряд ли появится доступ, так что о нем можно пока забыть. И, кстати, я и не предполагала, что он до такой степени неповоротливый.

Так что, если воспроизводимость в общих чертах есть, то, наверное, можно этим и ограничиться.

Да, мне не нужны точные совпадения по энергиям, важен общий вид зависимостей.
Мне очень понравилась предложение sanya1024 взять метод CCSD(T) в качестве эталонного, для этого либо считаю в GAMESS-US, либо прошу у Грановского FF 8. Но это, видимо, будет позже.

Сейчас провожу все расчета в базисе 6-31G(d)+. Дальше будет видно.

По ацетону в газовой фазе результаты уже выкладывала. Но вот, что получилось дополнительно: как и предполагала sanya1024

в долине реагентов есть минимум, найденный авторами, а рядом, отделенный небольшим барьером, лежит другой минимум, поглубже.

И этот минимум очень мал или его вообще нет. При уменьшении шага (stride) IRC c 0.5 до 0.3 и ужесточении условий сходимости (orttol=0.00001) получается следующая картина:

IRC-6.gif

Длина связи 4,02 это как раз геометрия с отъехавшей боргидридной группой.
Подобное же наблюдается и в случаях с добавлением растворителя, как для ацетона, так и для 4-метилциклогексанона. Результаты выложу в следующих постах.
Изначально я думала, для того, чтобы проскочить этот небольшой барьерчик, нужно увеличивать stride и уменьшать orttol, но я ошиблась.

А $FORCE NVIB=2 PURIFY=.t. $END, действительно очень помогает.

[Shadin]

Восстановление 4-метилциклогексанона, аксиальная атака (учет растворителя). Здесь все получилось сразу (STRIDE=0.5 и OPTTOL=0.0001), из ПС к реагентам:

IRC.gif

Именно это и заставило меня сделать расчеты, результаты которых приведены в предыдущем сообщении.

[Shadin]

Восстановление 4-метилциклогексанона (экваториальная атака) в растворителе еще считаю.
По ацетону могу сказать, что на пути от реагентов к ПС имеется еще одно ПС (ПС-1) с r(C-H)=1.49 ангстрем:

ацетон в р-ре.PNG

Зависимость привожу от длины связи, так как IRC из ПС останавливается на структуре с r=1.38, дальше был поиск ПС-1 и IRC уже из него. Здесь суммарная кривая.

[o-oxhem]

Значит, там есть тока бифуркации valley ridge inflection point (VRI). Её можно найти. Т.е. реально ПМЭР отклоняется от кривой IRC, в этом случае IRC яляется плохой моделью пути реакции. Тема обсуждалась на форуме. Да, и вообще надо чтить работы, например, Базилевского, особенно, когда им удавалось опередить ученых США, Европы и пр. и всячески это подчеркивать.

PS: Надеюсь, что Вы проверили, что точка, именованная как ПС-1, действительно переходное состояние. Можно запустить IRC, тогда Вы узнаете реальные долины реагентов (или продуктов, смотря откуда смотреть), к которым ведет ПМЭР.

[sanya1024]

Добрый день всем! я уезжала на несколько дней, возвращаюсь -- а тут такая бурная активность

Добавлю маленькое техническое замечание: Coupled Clusters работают в GAMESS, а в FireFly они если и есть, то в зачаточном состоянии.
Кстати, добыла бесплатную и параллельную программу Psi4 (http://www.psicode.org/). Пока не тестировала на реальных системах, но выглядит привлекательно