Выбор метода для расчета энергии МО

[dycost]

Пытаюсь подобрать боле-менее экономичный метод для расчета энергий LUMO/HOMO, нужны для оценки возможности взаимодействия по теории граничных орбиталей.
Считаю относительно простые молекулы (H₂O, CH₄,H₂CO,CH₃NH₂).
Методы - DFT, как чистые так гибриды и даблгибриды и Хартри-Фок.
Сравниваю с потенциалами ионизации и сродству к электронам (если известны), ну или грубо по энергии LUMO нуклеофил/электрофил.
Сводная таблица для воды в прикрепленном файле.
Вопросов основных несколько.
1) Существуют ли функционалы из распространенных (считаю в Firefly) "заточенные" под энергию граничных орбиталей. Гугление дало статью Orbital energies and negative electron affinities from density functional theory: Insight from the integer discontinuity The Journal of Chemical Physics 129, 044110 (2008); Andrew M. Teale, Frank De Proft, and David J. Tozer, но сих функционалов в Firefly нет, и цели немного другие, у авторов было описать наблюдаемые значения сродства к электрону и первого потенциала ионизации. У наиболее распространенных B3LYP, PBE и проч. отклонения существенны, даблгибриды ближе, но совпадения с экспериментом скорее случайны.
2) С нуклеофилами вроде все получается, голый Хартри-Фок и базисы с диффузиоными функциями (aug-cc-pVDZ) дают приемлемые результаты. С электрофилами те же методы не срабатывают, энергия LUMO > 0 (формальдегид aug-cc-pVQZ HOMU/LUMO -0.4425/0.0235), это ограничения ХФ? Либо я что-то делаю не так?
И на том же формальдегиде что-то творится с формой облаков LUMO. Файлы прикреплены. Вкратце -- без диффузионных функций вроде все норм, как в учебниках, добавляем диффузионные функции получаем совсем иной вид + коэффициенты разложения перед АО > 1. С чем данное явление может быть связано? Меня смущает, то что коэффициенты разложения по АО больше единицы.
Входные и выходные файлы расчета могу приложить, если нужно.

[Гесс]

Не совсем по теме, но равенство энергии ХОМО энергии ионизации, а ЛУМО - сродству к электрону - это приближение, которое выполняется далеко не всегда. С другой строны HOMO-LUMO gap зачастую коррелирует с электрохимическим гэпом заметно лучше чем расчетные IP-EA. Там как то отыгрывается взаимокомпенсация ошибок, но деталей не знаю, да и к поставленному вопросу это не очень относится.
Касательно вида орбиталей - они настолько непохожи что у меня подозрение что это разные орбитали. Попробуйте посмотреть соответствия не только ЛУМО-ЛУМО разных методов а скаэем по 5 низших орбиталей в обоих базисах.

[dycost]

Попробуйте посмотреть соответствия не только ЛУМО-ЛУМО разных методов а скаэем по 5 низших орбиталей в обоих базисах.

Да, про это я не подумал.
В базисе без диффузионных функций ничего похожего на "бублик" нет.
А вот в базисе с диффузионными нечто похожее попадается.

В этом похоже разобрался.
По цифрам в LUMO/aug-cc-pVQZ дают существенный вклад 6s АО углерода и кислорода, которых просто нет в базисе cc-pVQZ, и похожая по остальным вкладам от АО орбиталь LUMO-2/cc-pVQZ есть LUMO-1/aug-cc-pVQZ, но выглядит немного по другому. Спасибо за наводку)

Остается понять почему энергия LUMO электрофила положительная остается, как ни увеличивай базис...

[Гесс]

1) Существуют ли функционалы из распространенных (считаю в Firefly) "заточенные" под энергию граничных орбиталей.

Короткая версия рожденная после написания длинной: вам надо уходить в full-long-range exchange функционалы которых в незаточенной под ДФТ Firefly нету.

Длинная версия:

Конкретно для энергии LUMO попробуйте следующие функционалы:
M06HF (он будет провираться в HOMO, учтите что M06, M06L, M06-2X это вообще далеко не то же самое).
wB97 и его друга wB97X
M11 (малопопулярное добро)
LC-BLYP и LC-wPBE (коль скоро файерфлай наследует гамесс, а в гамессе LC есть, то есть надежда)

Дальнейшие улучшения энергий орбиталей в рамках ДФТ будут включать в себя IP-tuning (искать работы Thomas Körzdörfer где то из района 10.1021/acs.jctc.5b00873 и на пару лет вокруг), и замешивание туда функций Грина (G0W0) (НЁХ)
Наконец пока вы играете в молекулы такого размера то можно посмотреть что получается из CC2 и ADC (algebraic diagrammatic construction) (еще одна НЁХ). ADC(3) должен работать очень хорошо для предсказаний IP, EA и gap (что там с орбиталями я не очень в курсе), но бесплатный софт имеет максимум ADC(2). Ни того ни другого в файерфлае нет, возможно есть в PSI4, но это серьезно перекраивание задачи. и будет совершенно не экономично.

Если ничего из этого нет в файерфлае или то что есть не работает - попробуйте бесплатную ORCA.

UPD - я посмотрел http://classic.chem.msu.su/gran/gamess/ ... _DFT-D.pdf - кажется Грановский решил с range-separated функционалами не заморачиваться и Трулара тоже проигнорил. Впрочем Firefly всетаки прога под другие цели, так что это нормально. Таки в орку. Ну или в PSI4.

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

[mizuchi]

Firefly выдает какой-то ад. Посчитайте с помощью ORCA функционалом B2P-LYP в базисе def2-TZVP.

[dycost]

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

Задача -- предсказать хемосорбцию на поверхности. Комплексы достаточно большие (>=100 атомов сорбат и столько же примерно сорбент), поэтому моделирование "в лоб" задача неподъемная пока. А вот отдельно посчитать сорбат и сорбент вполне можно. Пытаюсь через различные индексы реакционной способности оценить возможность и предполагаемую силу взаимодействия.
ИРС пока такие: энергия и локализация в пространстве граничных орбиталей, заряды на атомах (по Бейдеру) и карты молекулярного электростатического потенциала.
Подход примерно такой:
- считаем энергию орбиталей, смотрим "мягкие" или "жесткие" вещества
- если мягкие, смотрим возможно ли взаимодействие путем перекрывания орбиталей, по геометрии и энергии
- если жесткие, смотрим возможно ли электростатическое взаимодействие, по геометрии и зарядам
- ну и до кучи проверяем электрофильность/нуклеофильность отдельных участков молекул по картам МЭП, что тоже идет в оценку электростатического взаимодействия.
Немного вот тут писал, на примере фталоцианина меди: https://butlerov.com/files/reports/2005 ... 5-7-22.pdf

IP и EA предсказывать не нужно, я просто их использовал, чтобы хоть как-то оценить адекватность расчета. IP кстати примерно совпадает с HOMO (HOMO>IP), а вот LUMO оценить так и не удается, только по знаку.
Про предсказания IP и EA есть неплохая работа

1.2961035.pdf

, оттуда, кстати, брал некоторые значения энергий орбиталей для сравнения.

По поводу ORCA и PSI4 буду смотреть, линь они вроде поддерживают оба.

Посчитайте с помощью ORCA функционалом B2P-LYP в базисе def2-TZVP

B2PLYP я пробовал в Firefly, в первом сообщении в таблице сравнительной он есть, результаты так-себе.
def2-TZVP насколько я понял без диффузионных функций?
К сожалению, как показали расчеты, LUMO завышенным получается без диффузионных функций, на всех функционалах.

[Гесс]

Да весь квантовый софт на линь. Ну ок, орка идет на всем.

В зависимости от наличных ресурсов и использования ri-ускорялок 200 атомов это вполне подьемная задача чтобы считать без сложных моделей, надежность которых мне не очевидна.
Хемисорбция это образование химической связи. Во первых мне сомнительно что прочность связи удастся оценить из энергий орбиталей, во вторых дисперсионка будет все равно практически полностью проигнорена, а в комплексах она может быть очень большой.
На мой взгляд проблему надо решать в лоб
Максимальный вариант: DLPNO-CCSD//DFT+D3 с триплзета базисами в обоих случаях. По возможности конечно TightPNO.
Приемлимый вариант: RI-PBE, возможно с D3, базис по обстановке. Попробуйте скормить свой комплекс в Priroda, она быстра.
Если все ну совсем плохо: DFTB, GFN-XTB2 или ANI-1 если система ограничена HCNO элементами. Дисперсионка скорее всего будет некорректной, но лучше чем ее полное игнорирование в приближении независимых мономеров.

Посчитайте с помощью ORCA функционалом B2P-LYP в базисе def2-TZVP

B2PLYP я пробовал в Firefly

Насколько я понимаю mizuchi подозревает неправильность орбиталей выдаваемых Firefly. Я с фаейрфлаем не работал, так что ничего сказать не могу.

[Гесс]

Немного вот тут писал, на примере фталоцианина меди: https://butlerov.com/files/reports/2005 ... 5-7-22.pdf

Гхм. Индексы Фукуи. Я уже не думал что их ктото еще использует.
Судя по благодарностям у них был доступ к Ломоносому, и при этом все что они сделали в отношении машиновремени это была B3LYP/6-31G(d,p) оптимизация молекулы на полсотни атомов с офигенной симметрией. Даже в Firefly это должна быть задача подьемная на неигровом ноуте.

Впрочем еще раз: основная функция Firefly это мультиреференсные расчеты где она великолепна. DFT в ней - побочные игрушки, причем исходя из мануала - не особо впечатляющие. При этом из всего перечисленного выше софта (ORCA, Priroda, PSI4) у Firefly самый тяжелый формат инпута. Так что переход в более общие пакеты будет для вас не только в скорость, но, вероятно, и в user-friendliness.

[dycost]

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

Дисперсионная энергия это имеются ввиду силы Лондона?
Вещества в основном полярные рассматриваться будут, а она вроде как много ниже ориентационной и индукционной энергии. Но на всякий случай попробую анализ Морокумы для какой-нибудь из систем.
Насчет оценки прочности, нужна не точная энергия связи, а скорее ряд по прочности, имеем поверхность и N шт. молекул, которые на ней сорбируются, и отсортировать по предполагаемой прочности связи.

Попробуйте скормить свой комплекс в Priroda, она быстра

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

Буду пробовать ORCA. Скачал, как время будет начну курить маны.

[Гесс]

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

Но с другой стороны в периодике дешево тоже не получится. 100-атомная сорбирующаяся молекула потребует порядка 1нм*1нм поверхности если не больше, плюс как минимум 2 нм вакуума от поверхности, И уже с тремя слоями рутила это на ДФТ будет вылетать в очень нехреновую копеечку. Впрочем коль скоро подложка это оксид, то возможно ее сможет описать периодическая полуэмпирика, типа PM6 в Empire (бесплатно на одну ноду для академических пользователей). Там заодно и всевозможные локальные свойства типа MEP, EA, IP, в чем то близко к индексам Фукуи только с покером и девушками. Качество может подхрамывать, но цена безусловно приятная, так что если системы родственные то системная ошибка будет одинакова и тренд построить можно будет.

[mizuchi]

Соглашаюсь с Гессом по всем пунктам. Если бы задача заключалась только в расчете свойств простых молекул, то тяжелая ТФП, типа предложенных двойных гибридов, или методы, явно учитывающие корреляцию, имели бы смысл, но простая молекула плюс поверхность -- совсем другая история. Чтобы вытянуть стоатомный кластер через GTO, нужен суперкомпьютер. Если такового нет, начинайте курить маны чего-нибудь плосковолнового.

[dycost]

Для неорганики -- кластеры ~100 атомов, для органики 2-4 молекулы, те же 100+ атомов.
Насчет полуэмпирики... я пробовал MOPAC PM6-7, на кластере, без периодики, результаты для неорганики не очень и по геометрии и по энергии, та же сорбция воды на (110) поверхность вышла ~15 кДж/моль, эскперименты 60-100, различные DFT с LANL2DZ на металле и попловские двухэкспонециальные на кислороде с учетом BSSE около 65-90.

На базисе я пока остановился aug-cc-pVDZ, при увеличении до TZ и QZ энергия орбиталей мало меняется. Их я брал на пробу, для малых молекул.

С рутилом мне повезло, нужная орбиталь (LUMO, он явный элекрофил) оказалась на поверхности (110).
Пока буду читать и думать, инфы много, спасибо)

[Гесс]

https://www.researchgate.net/publicatio ... ns_Studies

Я статью почитал по диагонали, но там фигурируют энергии орбиталей на B3LYP в гауссиане, индексы Фукуи, какие то выводы из всего этого, правда кажется только для сорбата, поверхность там возникает только в монтекарло в COMPASSовском силовом поле, в MaterialStudio, которое 1) не бесплатное, 2) уже 2 года не Акцерлис, а Биовиа, 3) но есть на торрентах (не будем тыкать пальцами )

[dycost]

Да, там только для сорбата, для сорбента монте-карло + DFT из литературных данных (вроде только работа выхода электрона для железа). B3LYP/6-31G(d,p) я считал энергии орбиталей, выходит далеко от реальности, а вот фукуи может и прокатить. Разность электронной плотности, если считать одним методом для нейтральной молекулы и ионов -- может что-то и выйдет дельное.
За наводку на MaterialStudio отдельное спасибо. Молекулярную механику я пробовал на том же рутиле в поле Amber, завышало энергия до 300-400 кДж/моль, но в этой статье вроде цифры адекватные выходят. Попробуем)

P.S.
Кстати Фукуи...
Я их в Gabedit считал, и там не требуется ЭП катиона и аниона, только орбитали для нейтральной молекулы, я так понял, что алгоритм просто убирает/добавляет электрон HOMO/LUMO, без учета перераспределения электронной плотности при переходе к ионам. Насколько точным будет данное приближение? особенно для молекул, где HOMO находиться на ароматическом кольце.

[Гесс]

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

Посмотрите статьи Timothy Clark где он обсуждает локальные свойства типа EAL, IPL, MEP(L)

Молекулярную механику я пробовал на том же рутиле в поле Amber, завышало энергия до 300-400 кДж/моль

GAFF неплох, но там нет титана.