Моделирование никелевого катализатора

[Vit Nhoc]

У нас есть много экспериментальных данных по никелевому (а также медному) скелетному катализатору: его активность с различными адсорбатами (константы скорости реакций на его поверхности), и ИК-спектры адсорбатов над этим катализатором. Мы хотим посчитать всё это квантовохимически. Я полагаю, надо для начала взять молекулы Ni2, Ni3, Ni4 и т.д., и посмотреть как эти молекулы взаимодействуют с нашими адсорбатами (CO, CO2, C6H6 и пр.).
При попытке рассчитать молекулу Ni2 обнаружились большие проблемы. Простой расчёт методом B3LYP/6-311G(D,P) в Gaussian даёт нестабильное состояние (это обнаруживается расчётом STABLE). Можно получить результат чуть получше, если сначала запустить STABLE(OPT), а на следующем шаге обычный расчёт OPT FREQ с опцией GUESS(READ) GEOM(ALLCHECK). Я пробовал мультиплетность 1, 3, 5, энергия получилась наименьшая с 3. Но всё равно, боюсь, этот расчёт далёк от реальности: длина связи Ni-Ni получилась 2.240 A, а в никелевом кристалле она вроде равна 2.65 A. Если заменить функционал B3LYP на PBEPBE, расстояние станет равно 2.239 A.
Я попробовал расчёт TD(NSTATES=12, ROOT=0), но может быть так и не разобрался как правильно его запускать: расчёт оптимизации с этим ключевым словом дал энергию на 28 кДж/моль выше (чем расчёт с Guess(read)), а длину связи 2.243 A. При этом такой расчёт занял в 15 раз больше времени. Может быть, в строке Total Energy, E(TD-HF/TD-KS) = приводится энергия не самого низколежащего состояния, а того которое нашёл первоначальный расчёт? Я даже пробовал запустить TD(ROOT=-2), но такой расчёт обрывается.
Расчёт с опцией SCF(QC) дал энергию ещё выше, чем обычный SCF.
Вообще, я плохо разбираюсь в расчётах таких систем, где нужны многоконфигурационные подходы. Я даже не знаю, TDDFT это одноконфигурационный или многоконфигурационный расчёт?
Мои вопросы:
1) Можно ли для нашей задачи обойтись методом TDDFT, или даже простым STABLE(OPT), если не нужны сильно точные цифры? Предполагается что мы будем строить корреляции, нам наверно достаточно полуколичественного результата. Или тут не обойтись без методов вроде MCSCF?
2) Почему расчёт TD дал энергию выше, чем простой SCF после STABLE(OPT)? Может быть, я неправильно его провёл?

[Гесс]

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

Я не возьмусь судить насколько подход через кластеры будет близок к желаемым данным, равно как и не уверен что учет волновой нестабильности даст чтото очень полезное. Я вообще не специалист в таких расчетах но предположительно раз вам надо "полуколичественно" можно попробовать остаться в рамках ДФТ и не лезть в мультиконфигурационные методы.
Возможно чтото полезное можно получить используя BrokenSymmetry, то есть указывая в инпуте guess=mix

[VTur]

Не знаю, на чем и чем Вы считаете, но B3LYP Вы взяли зря.
Попробуйте следующие функционалы: LDA, TPSS, TPSSH, revTPSS, HTBS, OPBE, PBEsol. Начните с LDA, потом OPBE, не понравится, используйте остальные, но они дороже.

[Vit Nhoc]

Я пробую обсчитать молекулу Ni3, постоянно не сходится scf. Иногда помогает SCF(QC), но он считает намного дольше. Насколько правилен такой подход - сначала считать обычным SCF, а если не сходится, запускать SCF(QC)?

Возможно чтото полезное можно получить используя BrokenSymmetry, то есть указывая в инпуте guess=mix

Я попробовал указывать NOSYMM, помогает. А какая разница между guess=mix и nosymm?

[Гесс]

Я пробую обсчитать молекулу Ni3, постоянно не сходится scf. Иногда помогает SCF(QC), но он считает намного дольше. Насколько правилен такой подход - сначала считать обычным SCF, а если не сходится, запускать SCF(QC)?

Использовать scf=xqc который первую половину расчета пытается свестись пообычному и если нет - то вторую половину добивает на qc.
Если это не поможет - можно пробовать повышать число шагов, но обычно если расчет неможет свестись даже через qc - все настолько плохо, что лечить надо уже не сходимость, а метод.

Я попробовал указывать NOSYMM, помогает. А какая разница между guess=mix и nosymm?

guess=mix перемешивает орбитали и "всреднем по больнице" улучает их расположение.
nosymm же делает расчет молекулы в C1 симметрии независимо от того какую симметрию вы имеете.
Гауссиан не повышает симметрию относительно исходной, поэтому если исходная структура не симметрична (Small random displacement) - Nosymm не должен играть никакой роли.
Вцелом это два несвязанных кейворда.
Кстати возможно стоит попробовать считать синглетные состояния Unrestricted или Restricted open-shell (RO). По второму ничего не скажу, анрестриктед для формально closed-shellнутых систем иногда очень славно. Иногда дурит, хотя формально должен давать тот же результат что и Restricted.
Unrestricted включается для систем с мультиплетностью больше 1 автоматически, а для синглетов - просто указанием U перед функционалом:

#P UPBE/aug-cc-pvdz scf=xqc Opt Freq Pop=None

Не расписывание заселенностей экономит время и много унылого места в аутпуте.

[Yurii]

На металлоорганику заточен Schrodinger: в гауссиане можно и пролететь.

P.S.
А вообще все они халтурщики: например, если иметь исходный код гауссиана, то на ха́суэллах его можно сделать раз в 5 быстрей, чем их продажный исполняемый код.