Отрицат. энергия возбужд. в B3LYP. Решение в CASSCF XMCQDPT

[sanya1024]

Посмотрела "неудачные" выдачи. Там орбитали получились откровенно плохими. Вместо пары-тройки доминирующих конфигураций (т.е. с коэф. хотя бы ~0.3) -- много какой-то мелкой шелупони с коэффициентами ~0.1-0.15. Орбитали по виду тоже плохие: среди pi орбиталей затесались 2 штуки sigma (связывающая и разрыхляющая). И боюсь, что насильственное запихивание вместо них в акт. пространство какой-либо пары pi орбиталей (путем той самой перестановки, о к-рой я рассказывала) не поможет. Хорошо бы взглянуть на орбитали из "удачного" расчета (я обратила внимание, что это GUGA, а не ALDET, но по идее при достаточно хороших критериях точности результат отличаться не должен).

Кстати, о п. 4. Я очень люблю рисовать такие картинки, откуда возникает ICT состояние. Пусть у нас есть состояние S0, первое возб. состояние (скажем, pi->pi* HOMO->LUMO переход) со своим минимумом, называемым LE (Locally Excited), и еще какое-то состояние (например, HOMO-2->LUMO), вероятно, с переносом заряда (ICT, Intramolecular Charge Transfer). У него свой минимум.

No_ICT.gif

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

ICT.gif

Т.е., в окрестности квазипересечения состояния S1 и Sn (это необязательно S2) будут смешанного характера, а после него либо S1 станет HOMO-2->LUMO, либо HOMO и HOMO-2 поменяются местами (и тогда S1 формально останется HOMO->LUMO, но HOMO уже будет другой). И на кривой S1 появится второй минимум -- ICT. Понятно, что акт. пространство для такого процесса должно включать все необходимые орбитали. А если между S1 и Sn есть еще какие-то состояния, то и их тоже. Т.е., простая на первый взгляд задача может вдруг оказаться не по зубам никакому кластеру, когда разрастется CAS. А это может наступить неожиданно быстро.

И еще одно замечание. Если какой-то расчетный метод (например, CIS) неправильно воспроизводит относительное положение кривых S0 и Sn, мы можем вообще не получить никакого квазипересечения и, соответственно, не увидеть ICT. А в DFT вообще, варьируя вес ХФ обмена в функционале, можно получать любые картинки, от первой (без квазипересечения) до второй. Теперь понятно, почему DFT и CIS в таких вопросах доверия нет. Более того, даже CASSCF может врать в отн. положении кривых, и именно поэтому нужно либо MRCI (над CAS-ом), либо MS-PT.

[sanya1024]

А Грановский -- ну что сказать, гений, да

[lmihael]

Т.е., в окрестности квазипересечения состояния S1 и Sn (это необязательно S2) будут смешанного характера, а после него либо S1 станет HOMO-2->LUMO, либо HOMO и HOMO-2 поменяются местами (и тогда S1 формально останется HOMO->LUMO, но HOMO уже будет другой). И на кривой S1 появится второй минимум -- ICT. Понятно, что акт. пространство для такого процесса должно включать все необходимые орбитали. А если между S1 и Sn есть еще какие-то состояния, то и их тоже. Т.е., простая на первый взгляд задача может вдруг оказаться не по зубам никакому кластеру, когда разрастется CAS. А это может наступить неожиданно быстро.

Столкнулся с проблемой в точке с углом 90 градусов. При малом вкладе S0 вообще не хочет сходится. Попробовал повторить GUGA SA с весами 1,1 Вначале все шло хорошо. Потом пошли выбросы. Стал уменьшать DXMAX.

p1_pt90mc_opt.rar

И вот что заметил S1 и S0 периодически меняются местами!!! Этим и обусловлены выбросы. Т.е. старый S1 оптимизировался так, что стал ниже по энергии, чем старый S0. Вот пример оптимизации и сингл-поинт расчетов некоторых точек на кривой оптимизации.

p1_pt90mc_p5.rar

p1_pt90mc_p6.rar

p1_pt90mc_p3.rar

Как я понимаю этой же проблемой обусловлена несходимость SS-расчета и SA-расчета с малым вкладом S0 - непосредственно при оптимизации орбиталей S1 становится ниже S0. Обратите внимание на характерную несходимость в течении 40 итераций в точке предшествующей минимуму (в расчете с оптимизацией геометрии).

К сожалению, не понял, как контролировать заселенность орбиталей в процессе оптимизации. Таблица вида:

STATE # 1 ENERGY = -1214.197287522

CSF COEF OCCUPANCY (IGNORING CORE)
--- ---- --------- --------- -----
1 -0.632553 22220000
31 -0.066604 21120011
60 0.093743 22020020
124 0.055333 11220101
150 0.051205 12120110
325 0.071918 02220200
491 0.064195 22211000

выдается только в начале расчета .
Попробую оптимизировать в соответственных точках S0, но боюсь орбитали "убегут" - сталкивался с таким при CIS расчете. Как мне кажется оптимальным вариантом было бы оптимизировать S0 геометрию из точки 5 выкинутым или замороженным 22220000 состоянием. Можно ли как-нибудь это сделать в Firefly?

Не понятно как в таком случае определять щель между состояниями. И каков вообще физический смысл смысл получаемого таким образом состояния S1 (бывшего S0, с основным вкладом 22220000)? Ведь если мы будем оптиизировать его геометрию (даже с замороженным углом) она все равно опустится ниже состояния с основным вкладом 22211000? Эта разность и будет реальная щель?

[lmihael]

Может быть просто заморозить бывшую HOMO-орбиталь (72-ую) когда на ней минимум электронной плотности?

[sanya1024]

И вот что заметил S1 и S0 периодически меняются местами!!! Этим и обусловлены выбросы. Т.е. старый S1 оптимизировался так, что стал ниже по энергии, чем старый S0.
Как я понимаю этой же проблемой обусловлена несходимость SS-расчета и SA-расчета с малым вкладом S0 - непосредственно при оптимизации орбиталей S1 становится ниже S0. Обратите внимание на характерную несходимость в течении 40 итераций в точке предшествующей минимуму (в расчете с оптимизацией геометрии).

Надо включать state-tracking $mcscf ntrack=сколько_состояний_отслеживать istate=какое_состояние_оптимизируете $end

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

К сожалению, не понял, как контролировать заселенность орбиталей в процессе оптимизации. Таблица вида:

STATE # 1 ENERGY = -1214.197287522

CSF COEF OCCUPANCY (IGNORING CORE)
--- ---- --------- --------- -----
1 -0.632553 22220000
31 -0.066604 21120011
60 0.093743 22020020
124 0.055333 11220101
150 0.051205 12120110
325 0.071918 02220200
491 0.064195 22211000

выдается только в начале расчета .

Тут я тоже не знаю, можно спросить Грановского на форуме FF, есть ли возможность выдавать такую табличку на каждой геом. итерации.

Попробую оптимизировать в соответственных точках S0, но боюсь орбитали "убегут" - сталкивался с таким при CIS расчете. Как мне кажется оптимальным вариантом было бы оптимизировать S0 геометрию из точки 5 выкинутым или замороженным 22220000 состоянием. Можно ли как-нибудь это сделать в Firefly?

Попробуйте пока state-tracking, если не поможет, будем дальше думать.

А от "убегания" орбиталей хорошего средства нет. Расширение акт. пространства тут может даже навредить, а не помочь (хотя попробовать можно). Еще можно включить в усреднение больше состояний (если размер акт. пространства это позволяет). Например, в Вашем случае (вспомните картинку про ICT, ту, где есть псевдопересечение) можно было бы попробовать включить в усреднение S0, S1 и S2. До барьера в S1 доминирует одна конфигурация, а после -- другая. Причем я даже скорее имею в виду не последовательность двоечек, единичек и ноликов, а тип орбиталей, откуда и куда идет переход. И хорошо бы обе эти конфигурации учесть при движении вдоль координаты, а то неаккуратненько как-то Т.е., ставим nstate=3 (хочется надеяться, что там все-таки S2, а не какое-то высокое Sn), в усреднение с какими-то весами включаем все три состояния, и ntrack тоже ставим 3.

[lmihael]

А от "убегания" орбиталей хорошего средства нет. Расширение акт. пространства тут может даже навредить, а не помочь (хотя попробовать можно). Еще можно включить в усреднение больше состояний (если размер акт. пространства это позволяет). Например, в Вашем случае (вспомните картинку про ICT, ту, где есть псевдопересечение) можно было бы попробовать включить в усреднение S0, S1 и S2. До барьера в S1 доминирует одна конфигурация, а после -- другая. Причем я даже скорее имею в виду не последовательность двоечек, единичек и ноликов, а тип орбиталей, откуда и куда идет переход. И хорошо бы обе эти конфигурации учесть при движении вдоль координаты, а то неаккуратненько как-то Т.е., ставим nstate=3 (хочется надеяться, что там все-таки S2, а не какое-то высокое Sn), в усреднение с какими-то весами включаем все три состояния, и ntrack тоже ставим 3.

Попробую так сделать. Тут еще столкнулся с трудностью-неприятностью. Получил каким-то образом (помоему SA оптимизацией с ALDET и весами 1,1,1) геометрию и орбитали даюшие энергии S0 и S1 ~-1214.240 и -1214.160

p1_pt00mca_sa.rar

Потом беру эту геометрию, считаю орбитали с нуля HF далее CASSCF энергию и получаются энергии ~-1214.219 и -1214.141.

p1_pt00mc_ad.rar

Можно ли сказать что результат ~-1214.219 и -1214.141 заведомо неправильный, т.к. он выше по энергии? Или все-таки надо смотреть на орбитали? По орбиталям мне больше нравится последний результат.

[sanya1024]

А мне не очень нравятся оба варианта. Я смотрю на натуральные орбитали CASSCF. Пара орбиталей, задействованных в переходе (с заселенностями 1.3 и 0.7) выглядит в обоих расчетах почти одинаково. Это хорошо. Есть еще пара, к-рые можно как-то сопоставить. В выдаче *sa.* это орбитали с заселенностями 1.878 и 0.0977, а во втором это орбитали с заселенностями 1.8988 и 0.0858. Эти пары тоже похожи в обоих выдачах. А вот остальные орбитали не похожи совсем. В выдаче *sa.* есть еще пара pi орбиталей, локализованная на кольце C1-C2-C3-C4-C5-C6, и последняя пара -- sigma орбитали связи C3-C4 (это-то мне и не нравится). А в выдаче *ad*.* две оставшиеся пары орбиталей локализованы вообще на другом кольце. И это мне не нравится гораздо больше.

Вроде как расширять акт. пространство в этой задаче не очень хочется, но побороться за то, чтобы в акт. пространстве были только pi орбитали кольца C1-C2-C3-C4-C5-C6 можно. Мне почему-то кажется, что попадание в акт. пространство sigma орбитали связи C3-C4 -- случайное явление, но если в качестве старта для всех точек потенциальной кривой используются однотипные орбитали (полученные одним и тем же способом или орбитали от соседней точки), то эта случайно попавшая орбиталь будет кочевать из точки в точку. Так что можно попробовать считать орбитали из панча от *sa.*, поискать pi орбиталь нужного кольца под акт. пространством (вряд ли она убежала куда-то глубоко) и комплементарную ей над акт. пространством, ручками переставить (думаю, разберетесь) и запустить CAS с такого старта. Если все-таки эта сигма проявит упорство -- значит, судьба. Но можно будет еще попробовать вариант с включением в усреднение второго синглета, S2. Если и после этого упрямая сигма останется -- значит, точно судьба

[sanya1024]

Сейчас полезла искать орбитали на замену этой сигме -- фигушки! ничего подходящего в пределах 14 орбиталей в ту и другую сторону. То ли вправду судьба, то ли стоит попробовать включить в усреднение S2 (про триплеты временно забудем, делаем GUGA расчет или ALDET с pures=.t.)

[Ferom]

Можно ли сказать что результат ~-1214.219 и -1214.141 заведомо неправильный, т.к. он выше по энергии? Или все-таки надо смотреть на орбитали? По орбиталям мне больше нравится последний результат.

А идея Olion про синглетный бирадикал не вариант?

[sanya1024]

А идея Olion про синглетный бирадикал не вариант?

Ну, при повороте на 90 град. наш open-shell singlet (тот, где ..11..) вроде как и есть синглетный бирадикал: одна единичка соответствует орбитали, локализованной на кольце C1-C2-C3-C4-C5-C6, а вторая -- повернутой CHOH группе. А в плоской структуре такой четкой локализации нет.

[Ferom]

Столкнулся с проблемой в точке с углом 90 градусов. При малом вкладе S0 вообще не хочет сходится.

Недавно делал расчеты контактной пары ион радикалов в UHF с мультиплетностью 3. Сходимость была ужасная, получались то положительные, то отрицательные выбросы и так до бесконечности. Я так понял, что на одной итерации получалась конфигурация пары триплет + синглет, а на другой радикал + радикал. Проблему решил тем, что нашел для одной из геометрий вектор для конфигурации радикал + радикал и сохранил его, а затем при расчете других геометрий делал MOREAD.

lmihael, если у Вас в других точках (80, 70, …) сходимость нормальная, то попробуйте для проблемной точки 90 в качестве начальных МО использовать вектор из предыдущих удачных расчетов, может быть этот способ и в Вашем случае сработает.

[sanya1024]

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

[Himera]

Разговоры об орбиталях отправлены в старую тему

[lmihael]

А мне не очень нравятся оба варианта. Я смотрю на натуральные орбитали CASSCF. Пара орбиталей, задействованных в переходе (с заселенностями 1.3 и 0.7) выглядит в обоих расчетах почти одинаково. Это хорошо. Есть еще пара, к-рые можно как-то сопоставить. В выдаче *sa.* это орбитали с заселенностями 1.878 и 0.0977, а во втором это орбитали с заселенностями 1.8988 и 0.0858. Эти пары тоже похожи в обоих выдачах. А вот остальные орбитали не похожи совсем. В выдаче *sa.* есть еще пара pi орбиталей, локализованная на кольце C1-C2-C3-C4-C5-C6, и последняя пара -- sigma орбитали связи C3-C4 (это-то мне и не нравится). А в выдаче *ad*.* две оставшиеся пары орбиталей локализованы вообще на другом кольце. И это мне не нравится гораздо больше.

Вроде как расширять акт. пространство в этой задаче не очень хочется, но побороться за то, чтобы в акт. пространстве были только pi орбитали кольца C1-C2-C3-C4-C5-C6 можно. Мне почему-то кажется, что попадание в акт. пространство sigma орбитали связи C3-C4 -- случайное явление, но если в качестве старта для всех точек потенциальной кривой используются однотипные орбитали (полученные одним и тем же способом или орбитали от соседней точки), то эта случайно попавшая орбиталь будет кочевать из точки в точку. Так что можно попробовать считать орбитали из панча от *sa.*, поискать pi орбиталь нужного кольца под акт. пространством (вряд ли она убежала куда-то глубоко) и комплементарную ей над акт. пространством, ручками переставить (думаю, разберетесь) и запустить CAS с такого старта. Если все-таки эта сигма проявит упорство -- значит, судьба. Но можно будет еще попробовать вариант с включением в усреднение второго синглета, S2. Если и после этого упрямая сигма останется -- значит, точно судьба

Сейчас полезла искать орбитали на замену этой сигме -- фигушки! ничего подходящего в пределах 14 орбиталей в ту и другую сторону. То ли вправду судьба, то ли стоит попробовать включить в усреднение S2 (про триплеты временно забудем, делаем GUGA расчет или ALDET с pures=.t.)

Пересчитал несколько точек 0 30 90 получается все как в ad, т.е. минимум 2 пары орбиталей локализованы на тозильном кольце . Энергия S1 состояния как в ad - порядка -1214.140. Сигма орбитали нет. Пересчитал с усреднением по 4-ем состояниям - тоже самое с орбиталями и энергией. Видимо доверять стоит "ad" результату.

Столкнулся с проблемой в точке с углом 90 градусов. При малом вкладе S0 вообще не хочет сходится.

Недавно делал расчеты контактной пары ион радикалов в UHF с мультиплетностью 3. Сходимость была ужасная, получались то положительные, то отрицательные выбросы и так до бесконечности.

Мне кажется в Вашем случае как-раз сработает state-tracking. А у меня в любом случае орбитали убегают при оптимизации геометрии.

[lmihael]

1. Хотел еще спросить про расчет спектров. Правильно ли я понимаю, что вначале следует провести оптимизацию геометрии S0 методом CASSCF SS, а затем провести SA (например по 4-м синглетным состояниям) расчет энергий + XMCQDPT? Или геометрию S0 лучше иначе оптимизировать?

2. Попробовал я state-tracking расчет, с ним идет без вылетов, но орбитали убегают . Ну и в конце все-таки произошел мегавылет - состояния все-таки поменялись местами, я в форуме firefly отписал по этому поводу.

3. Взял я ту точку 90 где энергии S1 и S0 практически равны и рассчитал XMCQDPT и вроде бы то что было S1 здорово опустилось ниже старого S0.

p1_pt90mc.rar

p1_pt90xmc.rar

Есть ли какой-нибудь надежный способ идентифицировать состояния после XMCQDPT?

4. Попробовал прикинуть спектр с помощью XMCQDPT получилось похоже на эксперимент и CIS*0.7 и B3LYP с одним кардинальным отличием. И CIS и B3LYP давали два очень близко-лежащих состояния 321 и 319 нм с близкими силами осцилляторов. XMCQDPT в этой области дает только одно состояние!

5. Ну и наконец отнес письмо счастья в НИВЦ МГУ, так что надеюсь скоро мне разрешат доступ на кластер. И моя машинка отдохнет. Может на кластере взять базис помощнее?

[sanya1024]

1. Хотел еще спросить про расчет спектров. Правильно ли я понимаю, что вначале следует провести оптимизацию геометрии S0 методом CASSCF SS, а затем провести SA (например по 4-м синглетным состояниям) расчет энергий + XMCQDPT? Или геометрию S0 лучше иначе оптимизировать?

S0 надо оптимизировать так же, как S1. Я бы даже сделала так: если при оптимизации S1 веса были 0.1 и 1.0, то для оптимизации S0 поставила бы 1.0 и 0.1. Т.е., все должно быть сделано на одном уровне. Тут ситуация как раз неприятна тем, что S0 и S1 не только квазивырождены, но и меняются местами. Поэтому тем более с ними надо обращаться одинаково.

2. Попробовал я state-tracking расчет, с ним идет без вылетов, но орбитали убегают . Ну и в конце все-таки произошел мегавылет - состояния все-таки поменялись местами, я в форуме firefly отписал по этому поводу.

Пока там ответа нет, я еще персонально попрошу Грановского глянуть на этот вопрос (возможно, он еще не успел заглянуть на форум).

3. Взял я ту точку 90 где энергии S1 и S0 практически равны и рассчитал XMCQDPT и вроде бы то что было S1 здорово опустилось ниже старого S0.

p1_pt90mc.rar

p1_pt90xmc.rar

Есть ли какой-нибудь надежный способ идентифицировать состояния после XMCQDPT?

Открываем выдачу ChemCraft-ом и -- параллельно -- смотрим на текст. После строчки *** MC-XQDPT2 ENERGIES *** идет табличка с энергиями, столбец 2ND ORDER -- то, что нам нужно. Посчитали щель между S0 и S1, получилось 0.135 эВ. Это действительно маленькая щель, что весьма способствует безызлучательным переходам. Смотрим табличку после строчки EIGENVECTORS OF THE EFFECTIVE HAMILTONIAN. Там, где на диагонали этой таблички число, близкое к единице, значит, соотв. CAS-овское состояние после QDPT осталось на месте. Какие-то состояния могли смешаться (числа в районе 0.3--0.7). А CAS-овские состояния S0 и S1 поменялись местами ("почти единички" уехали с диагонали).

Теперь смотрим ChemCraft-ом орбитали (State 1 QDPT natural orbitals, State 2 QDPT natural orbitals). State 1 QDPT natural orbital 4 и State 1 QDPT natural orbital 5 имеют заселенности ок. 1, т.е., нижнее состояние -- конфигурации ...11... А State 2 QDPT natural orbitals все имеют заселенности ок. 2, т.е., второе состояние -- конфигурации ...20... Вот так. Все-таки поменялись они местами. Кстати, 3-е состояние (S2), к-рое мы упорно игнорировали, похоже по орбиталям на новый S0 (конфигурации ...11...). Может, правда, стоило его учитывать в CAS-е? Если посмотрите на картинку ICT.gif несколькими постами выше, то увидите, что где-то в середине потенциальной кривой оно подходит близко и должно давать существенный вклад в S1, поэтому неплохо бы присматривать за ним (пусть даже не оптимизируя) по всей кривой, т.е., включить в усреднение и откладывать на графике.

Теперь смотрим после строчки RADIATIVE TRANSITION MOMENT BETWEEN STATE # 1 AND STATE #
Из TRANSITION DIPOLE (столбец NORM) вытаскиваем дипольные моменты переходов, пересчитываем их в силы осцилляторов -- и спектр готов.

4. Попробовал прикинуть спектр с помощью XMCQDPT получилось похоже на эксперимент и CIS*0.7 и B3LYP с одним кардинальным отличием. И CIS и B3LYP давали два очень близко-лежащих состояния 321 и 319 нм с близкими силами осцилляторов. XMCQDPT в этой области дает только одно состояние!

Очень может быть. Просто замечательно! В QDPT мы учитываем динамическую корреляцию, ту, что уже включена в корреляционный функционал, но полностью отсутствует в CAS. Зато в CAS есть статическая корреляция (та самая, что работает в открытых оболочках и отсутствует в DFT). Таким образом, после QDPT мы более-менее учитываем всю корреляцию, и все состояния встают на место. Именно поэтому я рассматриваю QDPT как способ верификации результатов TDDFT.

Близкие переходы в реальном спектре скорее всего сольются в один, а в QDPT у этого состояния наверняка будет удвоенная сила осциллятора по сравнению с "расщепившимися" в TDDFT состояниями.

5. Ну и наконец отнес письмо счастья в НИВЦ МГУ, так что надеюсь скоро мне разрешат доступ на кластер. И моя машинка отдохнет. Может на кластере взять базис помощнее?

Можно попробовать Удачи! Но я бы на Вашем месте включила бы в рассмотрение S2.

[lmihael]

Но я бы на Вашем месте включила бы в рассмотрение S2.

Воспользовался Вашим советом и включил S2. Еще раз внимательно посмотрел на орбитали и озадачился вопросом: правильно ли я понимаю, что в данном примере имеет смысл 70-ую орбиталь заменить 67-ой, а 76-ую 77-ой?

[sanya1024]

Попробовать можно, но, боюсь, они "убегут". Еще показались подозрительными вот эти двукратно-возбужденные конфигурации:
703 0.270735 12211100
890 -0.210492 22201100
в S1 и вот эти
890 0.275483 22201100
893 0.358103 12211100
в S2. Великоват у них вес...

Ну, и теперь смотреть, как ведут себя все три состояния при откручивании CHOH группы. Причем бОльший интерес представляют не CAS-овские результаты, а QDPT-шные.

[sanya1024]

Как выбирать активное пространство в CASSCF: свежая статья

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

http://onlinelibrary.wiley.com/doi/10.1002/qua.23068/abstract

[krv-nn]

Как выбирать активное пространство в CASSCF: свежая статья

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

http://onlinelibrary.wiley.com/doi/10.1002/qua.23068/abstract

А Вы не могли бы выложить pdf этой статьи здесь, если у Вас есть доступ? Наш университет, к большому сожалению, не продлил подписку в т.ч. на Wiley. Буду очень благодарен!