оптимизация сложных случаев в Gaussian

[morua]

Прошу более опытных коллег поделиться подходами к оптимизации в Gaussian (G09) трудносходящихся систем. Уже довольно длительное время безрезультатно пытаюсь оптимизировать возбужденное состояние интересующей меня молекулы в вакууме:

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

%name.chk
#p opt=(ReadFC,tight) Int=UltraFine TD=(Read,NStates=6,Root=1) b3lyp/6-31g(d,p) Geom=Check Guess=Read)

Чтение инструкций, в частности на официальном сайте программы и тертьей главы книги “Exploring Chemistry With Electronic Structure Methods” помогает мало, поскольку применение опции Opt=CalcAll и для основного состояния довольно затратное, а для возбужденного — даже думать боюсь.

Из других советов пытался применить Opt=(...,MaxStep=10) и SCF=QC, но безрезультатно.

[Shorku]

Я не специалист по гауссиану и оптимизации возбужденных состояний, но все же, если подтянутся более продвинутые спецы, опишите, в чем проблема?
Это осцилляции? Или оптимизации выходит на пологое плато и ползет, ползет... Или расчет просто в какой-то момент падает? Вы пробовали менять координаты (внутренние на декатртовы и наоборот), или метод оптимизации, проверяли, не может ли там быть вырождения - т.е. не сползает ли система до пересечения? Разумно ли выглядят промедуточные геометрии, или кровькишкирас"№%ло?

[morua]

Спасибо большое за уточняющие вопросы!
(1) Это осцилляции. Причем амплитуда осцилляций зависит от начальной геометрии (ее я тоже пытался искусственно менять, поскольку небольшая деформация помогала при расчете более сложных аналогов проблемной молекулы (винилогов)) и параметра MaxStep.
(2) Координаты менять не пробовал.
(3) Что касается вырождения, то это наиболее тревожный для меня вопрос. Дело в том, что, согласно TDDFT расчету, первые два возбужденные состояния очень близки по энергии, хоть и различной природы (pi-pi* и n-pi*). И если в процессе оптимизации их относительное расположение поменяется, то ... дальше я не знаю, увы. Но, по-видимому, получится, что изменится целевое состояния, задаваемое параметром Root=1.

[Shorku]

Замена координат обычно помогает перейти от состояния "на первом шаге оптимизации получаем звездное небо" к "оно соптимизировалось", видимо, не Ваш случай. Насколько велики амплитуды осцилляций? Как изменяется ТDDFT спектр в процессе? Гауссиан выдает энергии состояний для промежуточных геометрий?

[morua]

Промежуточные энергии выдаются (пример ниже в коде). А вот изменение TDDFT спектра мне не нравятся. Дело в том, что в процессе оптимизации полоса интересующего меня длинноволнового перехода смещается всего на ~10 нм (от 402 до 412-413), тогда как для винилога исследуемого соединения (рассчет которого прошел без проблем) расчетный виниленовый сдвиг составляет аж 28 нм (445 и 473 нм).

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

 SCF Done:  E(RB3LYP) =  -1055.02071691     a.u. after    1 cycles (initial)
 SCF Done:  E(RB3LYP) =  -1055.01826443     a.u. after   28 cycles
 SCF Done:  E(RB3LYP) =  -1055.02069743     a.u. after   27 cycles
 SCF Done:  E(RB3LYP) =  -1055.02071688     a.u. after    2 cycles
 SCF Done:  E(RB3LYP) =  -1055.02071691     a.u. after    1 cycles
 SCF Done:  E(RB3LYP) =  -1055.01836165     a.u. after    5 cycles
 SCF Done:  E(RB3LYP) =  -1055.01804214     a.u. after    5 cycles
 SCF Done:  E(RB3LYP) =  -1055.01691185     a.u. after    5 cycles
... (some steps are omitted)
 Step number  27 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01508457     a.u. after    5 cycles
 Step number  28 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01718971     a.u. after    4 cycles
 Step number  29 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01503095     a.u. after    5 cycles
 Step number  30 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01508453     a.u. after    5 cycles
 Step number  31 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01718970     a.u. after    4 cycles
 Step number  32 out of a maximum of  258
 SCF Done:  E(RB3LYP) =  -1055.01503095     a.u. after    5 cycles

(the task was terminated (cf. the energies of 29th and 32nd steps)

А еще смущает то, что после выдода на “зацикливание” постоянно меняются местами два состояния:

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

 New state      1 was old state      5
 New state      5 was old state      1

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

[Shorku]

Я имел в виду тд-дфт энергии всех рассчитываемых состояний в каждой точке оптимизации. Можно было бы наглядно посмотреть, пересекается там что-нибудь, или нет. Хотя последнее Ваше замечание про смену порядковых номеров наводит на нехорошие мысли.
Т.е. гауссиан умеет отслеживать, какое состояние каким было? Это хорошо, я думал, такая штука только в FireFly есть. Если состояния все-таки пересекаются то сколько состояний ни считайте, проблема так и останется. Ежли оно не слишком большое и тяжелое, могу попробовать, например, оркой пооптимизировать.

[morua]

Спасибо большое за предложени насчет Orca. но молекула действительно немаленькая (C22NO2H19). Попытаю счастья еще Удалось расчитать эту же молекулу в поле растворителя (PCM), но там благодаря сольватации “перепутывающие” в приближении вакуума состояния значительно сильнее отличаются по энергии. Хочу теперь взять геометрию и гессиан и TD из той задачи, убрать сольватацию и надеяться на праздник.

[sanya1024]

Если это не тайна, не покажете молекулу? похоже, действительно, перепутываются близколежащие состояния. В FireFly для этого есть state-tracking, если хотите, я попробую посмотреть, как пойдет с трекингом. Несколько смущает, что глюки происходят с 1-м возб. состоянием (обычно оно хорошо отщепляется от остальных и ни с чем не путается, кроме основного, но тогда это система с (псевдо)пересечением уровней и, соответственно, прощай single-reference).

[morua]

Спасибо, sanya1024. Молекула в приложении (out гауссиана для основного состояния).
И, хоть это и не соотсествует теме: не подскажете, как “включать” state-tracking в FireFly.

И еще, позор на мою голову: не заметил вначале, что для этой молекулы первое возбужденное состояние является (n-pi*), а (pi-pi*) только второе.

Добавлю, что учет сольватации (PCM) привел к существенному понижению энергии (pi-pi*)-состояния (по сравнению с (n-pi*)), и оно стало S1. Может этим и объясняется то, что в этом случае оптимизация геометрии возбужденного состояния прошла успешно?

[morua]

Еще раз Большое спасибо всем, кто отозвался на мой вопрос. Но отчаянная попытка провести расчет исходя из геометрии (а также гессиана и силовых констант), полученной для сольватированной в этаноле (PCM) молекулы неожиданно увенчалась успехом. Конечно, еще нужно проверить результат, рассчитать колебательную задачу (freq) для S1-состояния, но свет уже забрезжил в конце тоннеля

[sanya1024]

Спасибо, sanya1024. Молекула в приложении (out гауссиана для основного состояния).
И, хоть это и не соотсествует теме: не подскажете, как “включать” state-tracking в FireFly.

Это проще всего: обычно указывают NSTATE=<кол-во состояний> ISTATE=<номер целевого состояния>. Так вот, если этот номер указать отрицательным (в CIS и TDDFT оптимизации), то включится трекинг.

И еще, позор на мою голову: не заметил вначале, что для этой молекулы первое возбужденное состояние является (n-pi*), а (pi-pi*) только второе.

А, это все объясняет... почти. Сейчас поставила считать и смотреть, как оно себя ведет, в самое ближайшее время расскажу.

Добавлю, что учет сольватации (PCM) привел к существенному понижению энергии (pi-pi*)-состояния (по сравнению с (n-pi*)), и оно стало S1. Может этим и объясняется то, что в этом случае оптимизация геометрии возбужденного состояния прошла успешно?

А это как раз очень типично для таких красителей. Именно так и получилось.

[sanya1024]

Итак, рассказываю. Простая оптимизация 1-го состояния (того самого, к-рое оказалось неправильным, n-pi*) с ISTATE=1 сошлась спокойно, без осцилляций (видимо, дефолтные алгоритмы оптимизации в FireFly и Гауссиане разные). Финальная энергия -1054.1310044274. Везде оптимизируемое состояние называется State 1, переход 86->88, n-pi*. Оптимизация тоже 1-го состояния, но с трекингом (ISTATE=-1) сошлась туда же: финальная энергия -1054.1310044274, оптимизируемое состояние State 1, переход 86->88, n-pi*. Т.е., для n-pi* состояния на ППЭ есть минимум.

А вот оптимизация 2-го состояния с трекингом (ISTATE=-2), того, что в исходной точке было pi-pi*, повела себя интересно: уже на 1-м (отсчет от 0-го) шагу оптимизации получилось следующее сообщение:
WARNING: STATE FLIPPING OCCURED.
FIREFLY HAS FOLLOWED STATE 2, NOW FOLLOWS STATE 1 WITH OVERLAP OF 0.978036
И дальше пошло оптимизироваться 1-е состояние, оставаясь при этом pi-pi* характера. За 1 час расчета прошло 10 итераций (32 процессора), градиент еще немаленький, и энергия этого состояния дошла до -1054.1388576, т.е. ниже, чем энергия n-pi* состояния в минимуме! Т.е., pi-pi* состояние где-то пересекается с n-pi* и становится 1-м. Оптимизация 1-го состояния без трекинга в FireFly до этого места не дошла: или остановилась раньше (порог сходимости по градиенту 1e-4), или просто не попала в эту область (помним про разные алгоритмы оптимизации). А Гауссиановская оптимизация, видимо, как раз попала в область пересечения и принялась скакать с терма на терм, не зная за какое состояние хвататься.

Когда дозреет оптимизация S2 c трекингом, выложу отдельно.

[sanya1024]

Ну вот, дозрела и вторая оптимизация, выкладываю. Обратите внимание: в первом файле стоит ISTATE=-2, но поскольку состояния поменялись местами и под конец 1-го файла оптимизировалось уже 1-е состояние, то во втором файле надо ставить ISTATE=-1: при рестарте программа кушает только то, что в инпуте, а в инпуте -- геометрия и орбитали с последней точки. И в этой геометрии с этими орбиталями наше целевое состояние -- 1-е.