Противоречие в описании вращения TDDFT и QDPT методами

[sanya1024]

Ага, попробуйте. В CASSCF "хороший" старт приводит к более низкой энергии. Вот только если у нас SA-CAS, то понижается энергия всей пачки усредняемых состояний, а каждое отдельное состояние -- как повезет. Поэтому второй критерий -- отчетливое доминирование одной конфигурации или двух (в окрестности (квази)пересечения) с примерно равными весами. Собс-но, никто не мешает в каждой точке брать в качестве старта те орбитали, к-рые дадут лучший результат. Иногда хорошо работают орбитали от соседней точки. Главное -- чтобы конфигурационный состав не дергался от точки к точке.
Еще мне показалось, что 13 состояний в усреднении -- это перебор. Мне показалось, что для начала можно попробовать 7 состояний, а дальше наращивать только если будут какие-то негладкости или резкие скачки. Зато проследить за конфигурационным составом всех усредняемых состояний на протяжении всей кривой.

[lmihael]

Попробовал родные TDDFT орбитали в каждой точке без перестановок. Активное пространство 8x8, усреднение по 7 синглетам.
Конфигурации скачут. Энергии тоже.

mcscf.jpg

QDPT.jpg

MCSCF_8x8_N7_wo_exch.rar

XMC.rar

[sanya1024]

Так я и думала. Это разрывы. До 45 град идет один набор конфигураций, а после -- другой. А в верхних состояниях вообще страх божий... Явно видны множественные (псевдо)пересечения ППЭ. А из-за усреднения многие состояния просто выпали из нашего поля зрения. Да, CASSCF -- это искусство... нет, адская кухня!
Я понимаю, что сейчас предлагаю шаманские рецепты, но... возьмите точку 50 град. и с ее орбиталей стартаните CASSCF для 45 и 60 град. Посмотрите, насколько изменится результат.
Сейчас внимательно посмотрю на выдачи, может, что-то более конкретное скажу.

[sanya1024]

Ну, все плохо уже с первой точки, 0 град. В состоянии S0 ведущая конфигурация 22202000 с весом 0.94 -- это никуда не годится.
Как-то по всему похоже, что акт. пространство (8,8) маловато. Добавьте 2 занятых орбитали и попробуйте снова (стартуя с TDDFT). Да, надеюсь, Вы геометрию не переоптимизируете каждый раз? нам пока надо найти подходящее акт. пространство и усреднение, а как найдем -- тогда можно и геометрию дооптимизировать в каждой точке.

[sanya1024]

Ой, я только сейчас заметила: Вы все время считываете norb=200 орбиталей. А надо -- 390. Для CASSCF это не очень важно, а в QDPT может оказаться критично.

[lmihael]

Так я и думала. Это разрывы. До 45 град идет один набор конфигураций, а после -- другой. А в верхних состояниях вообще страх божий... Явно видны множественные (псевдо)пересечения ППЭ. А из-за усреднения многие состояния просто выпали из нашего поля зрения. Да, CASSCF -- это искусство... нет, адская кухня!
Я понимаю, что сейчас предлагаю шаманские рецепты, но... возьмите точку 50 град. и с ее орбиталей стартаните CASSCF для 45 и 60 град. Посмотрите, насколько изменится результат.
Сейчас внимательно посмотрю на выдачи, может, что-то более конкретное скажу.

В MCSCF энергии выровнялись

MCSCF.jpg

MCSCF_45&60_from50.rar

В QDPT тоже.

QDPT.jpg

QDPT_45&60_from50.rar

Ну, все плохо уже с первой точки, 0 град. В состоянии S0 ведущая конфигурация 22202000 с весом 0.94 -- это никуда не годится.
Как-то по всему похоже, что акт. пространство (8,8) маловато. Добавьте 2 занятых орбитали и попробуйте снова (стартуя с TDDFT). Да, надеюсь, Вы геометрию не переоптимизируете каждый раз? нам пока надо найти подходящее акт. пространство и усреднение, а как найдем -- тогда можно и геометрию дооптимизировать в каждой точке.

Я уже пробовал ранее - все тоже самое. На всякий случай пересчитал:

r1_00.rar

Геометрии не пересчитываю, за исключением того случая когда добавлял точки 40,45,50,75 - пересчитал на всякий случай и 0,30,60,90 чтобы все были в равных условиях.

Ой, я только сейчас заметила: Вы все время считываете norb=200 орбиталей. А надо -- 390. Для CASSCF это не очень важно, а в QDPT может оказаться критично.

Честно говоря не думал, что такие высоколежащие орбитали играют какое-то значение. Или они имеют какой-то другой смысл?
Пересчитал с NORB=390 - результат тот же - максимальное изменение энергии 10^-12H

[sanya1024]

В MCSCF энергии выровнялись
В QDPT тоже.

У меня подозрение, что по каким-то причинам single-reference методы плохо описывают даже S0 в окрестности 0 град. Так что надо двигаться в том же направлении, от 50 град к 0. Получили точку 45 град -- берите ее орбитали и делайте точку 30 град. А с орбиталей от 30 -- точку 0. Будем надеяться, что при этом кривая будет более гладкой.

Я уже пробовал ранее - все тоже самое. На всякий случай пересчитал:

r1_00.rar

Я уже и сама взялась ковырять эту задачу. Убедилась, что какое акт. пространство и усреднение ни бери, в 0 град при старте с DFT-шных орбиталей упорно получается ведущая конфигурация ...202... в S0. Причем на вид орбитали какие надо. Если их поменять местами, то энергия S0 не меняется (как и должно быть), а ведущая конфигурация становится, конечно, ...220... При старте с ХФ орбиталей такого нет, ведущая конфигурация сразу получается ...220...
Я пошла с другого конца: от минимального акт. пространства (2,2) расширила до (8e,7o) -- переставляя орбитали, выкинула те, что локализованы только на бензольных кольцах. Усредняла только по S0 и S1 с намерением добавить другие состояния в усреднение потом. Смотрела пока только точки 0 и 90. Выкинутые орбитали не стремились вернуться в акт. пространство -- видимо, не больно-то им и хотелось Но энергия S1 при 90 град оказалась выше, чем в 0. Значит, надо ловить в усреднение то состояние, к-рое понижает энергию S1. При больших усреднениях мы его поймали случайно, но не опознали среди прочих.

А что было в CIS? меня, понятно, интересует не только S0 и S1, но и более высокие состояния. В частности, энергия какого состояния идет вниз при скручивании от 0 до 90. Именно это состояние надо постараться поймать в усреднение.

Геометрии не пересчитываю, за исключением того случая когда добавлял точки 40,45,50,75 - пересчитал на всякий случай и 0,30,60,90 чтобы все были в равных условиях.

Да, с геометрией -- все правильно. Получили точки -- теперь с ними экспериментируем.

Честно говоря не думал, что такие высоколежащие орбитали играют какое-то значение. Или они имеют какой-то другой смысл?
Пересчитал с NORB=390 - результат тот же - максимальное изменение энергии 10^-12H

Ну, в QDPT в принципе учитываются все возбуждения во внешнее пространство, так что лучше не рисковать.

И вот еще что: я же помню, что считала CAS/QDPT и для системы с 4-пиридиниевым гетероциклом, а не только с бензотиазолиевым. И вот в приаттаченной статье даже нашла упоминание об этом. А хорошенько покопавшись в своих файлах, нашла и сами расчеты... Считала я только при 0 и 90, был SA(4)-CASSCF(6e,7o) и X-QDPT. И получилось, что в скрученном состоянии в CAS щель между S0 и S1 0.5 эВ, а в QDPT 1.6 эВ, и абсолютная энергия S1 в скрученном состоянии ниже, чем в плоском. Тогда всю кривую я проходить поленилась, да и система была чуточку сложнее -- вместо NMe2 сидел азакраун-эфир. Так что мы добьемся нормальных кривых, я уверена.

[lmihael]

У меня подозрение, что по каким-то причинам single-reference методы плохо описывают даже S0 в окрестности 0 град. Так что надо двигаться в том же направлении, от 50 град к 0. Получили точку 45 град -- берите ее орбитали и делайте точку 30 град. А с орбиталей от 30 -- точку 0. Будем надеяться, что при этом кривая будет более гладкой.

Так на 50 градусах тоже не все ладно с конфигурацией:
STATE # 1 ENERGY = -764.758670902

CSF COEF OCCUPANCY (IGNORING CORE)
--- ---- --------- --------- -----
316 -0.115130 22200200
491 -0.114759 22211000
713 -0.067224 12201101
891 -0.050380 22111100
1275 0.944391 22202000
1276 -0.073473 22112000
1279 -0.117392 22022000
1310 -0.064415 02202002
1364 -0.077942 11202011

Я уже и сама взялась ковырять эту задачу. Убедилась, что какое акт. пространство и усреднение ни бери, в 0 град при старте с DFT-шных орбиталей упорно получается ведущая конфигурация ...202... в S0. Причем на вид орбитали какие надо. Если их поменять местами, то энергия S0 не меняется (как и должно быть), а ведущая конфигурация становится, конечно, ...220... При старте с ХФ орбиталей такого нет, ведущая конфигурация сразу получается ...220...
Я пошла с другого конца: от минимального акт. пространства (2,2) расширила до (8e,7o) -- переставляя орбитали, выкинула те, что локализованы только на бензольных кольцах. Усредняла только по S0 и S1 с намерением добавить другие состояния в усреднение потом. Смотрела пока только точки 0 и 90. Выкинутые орбитали не стремились вернуться в акт. пространство -- видимо, не больно-то им и хотелось Но энергия S1 при 90 град оказалась выше, чем в 0. Значит, надо ловить в усреднение то состояние, к-рое понижает энергию S1. При больших усреднениях мы его поймали случайно, но не опознали среди прочих.

Даже при больших усреднениях 90-градусов становится выгодно, только если орбитали при 0 будут "плохими" не выгодными.

А что было в CIS? меня, понятно, интересует не только S0 и S1, но и более высокие состояния. В частности, энергия какого состояния идет вниз при скручивании от 0 до 90. Именно это состояние надо постараться поймать в усреднение.

В CIS смотрел только S0 и S1 состояния. Сейчас посчитаю побольше.

И вот еще что: я же помню, что считала CAS/QDPT и для системы с 4-пиридиниевым гетероциклом, а не только с бензотиазолиевым. И вот в приаттаченной статье даже нашла упоминание об этом. А хорошенько покопавшись в своих файлах, нашла и сами расчеты... Считала я только при 0 и 90, был SA(4)-CASSCF(6e,7o) и X-QDPT. И получилось, что в скрученном состоянии в CAS щель между S0 и S1 0.5 эВ, а в QDPT 1.6 эВ, и абсолютная энергия S1 в скрученном состоянии ниже, чем в плоском. Тогда всю кривую я проходить поленилась, да и система была чуточку сложнее -- вместо NMe2 сидел азакраун-эфир. Так что мы добьемся нормальных кривых, я уверена.

Будем надеяться.

[lmihael]

Посчитал я CIS.

CIS_from MCSCF.rar

CIS.jpg

Если анализировать конфигурации то явно видно что S4 переходит в S3 в области 30^о, а новоявленный S3 уже в S2 в области 50^о

result.jpg

[lmihael]

Посчитал все точки последовательно с орбиталями из 50 (50-60-75-90 и 50-45-40-30-0).
Вот результат MCSCF:

MCSCF.jpg

MCSCF_from 50.rar

QDPT все равно не выгодно:

QDPT.jpg

QDPT.rar

[sanya1024]

Днем посмотрю файлы... пока по графику CASSCF видно, что кривая S1 какая-то слишком волнистая, что ли. Подозрительны мне лишние изгибы на ней. Мы добьем эту чертову систему, разберемся, что там не так но только днем, не сейчас...

[Alex Granovsky]

Всем доброй ночи!

А почему, собственно, эта связь (а также и другие связи на изначальном рисунке) должны так уж легко вращаться?
Если глянуть на дипольные моменты S0 и S1, при возбуждении происходит полная перезарядка, дипольный момент меняет направление. Как мне кажется, у S1 ведущая конфигурация та,
которая отвечет образованию хиноидных структур с заменой одинарных связей на двойные и наоборот. Это - в газе. Что в растворе, сказать намного сложнее.

[sanya1024]

Ну, из того, что я знаю про стириловые красители (а я знаю про них изрядно), как раз в основном состоянии они "хиноидные", а в возбужденном связи альтернируют. Поэтому такие стириловые красители слабо флуоресцируют (эксп. факт) -- флуоресценция тушится за счет вращения вокруг формально одинарных связей (кто бы мог подумать? )
С другой стороны, протонирование аминогруппы приводит к ровно таки обратному эффекту: в основном состоянии связи альтернируют, а в возбужденном краситель становится "хиноидным" -- с тенденцией к выравниванию длин связей, плоским и жестким как доска. И весьма интенсивно флуоресцирующим (не как родамин, конечно, но кв. выход уже измеримый получается). Наши экспериментаторы называют это катион-индуцированным разгоранием флуоресценции.
С третьей стороны, именно в протонированной форме центральная двойная связь (т.е., она была двойной в свободном красителе, а в протонированном стала какая-то "полуторная") становится более уязвимой, более подверженной фотохимическим реакциям в возбужденном состоянии -- транс-цис-фотоизомеризации или [2+2]-фотоциклоприсоединению. И эти процессы, конечно, иначе как в CASSCF/QDPT не сделаешь, TDDFT тут стопудово не проканает.
С полной перезарядкой при возбуждении согласна, есть такое дело. И это никак не противоречит вышесказанному.

[Alex Granovsky]

Беру выдачу для плоской структуры отсюда: download/file.php?id=116814

Для S0, если глянуть на заряды (малликеновские и левдиновкие, которые есть в выдаче, уж какие есть, при желании можно посчитать и другие) -

3 N 7.653074 -0.653074 7.008178 -0.008178
17 N 7.750844 -0.750844 7.142618 -0.142618

То есть азот на -N(CH3)2 конце заряжен отрицательнее

Для S1

3 N 7.698582 -0.698582 7.059317 -0.059317
17 N 7.701134 -0.701134 7.090229 -0.090229

и азоты почти выравниваются.

То есть все же, картина получается противоположной той, которую описывает sanya1024

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

[sanya1024]

Я, когда этим занималась (в основном на краунсодержащих красителях, но и на простых аналогах тоже), смотрела по длинам связей, а не по зарядам, поскольку с зарядами всегда есть к чему придраться. А структуры тут в файлах -- результаты релаксированных сканов для S1. Поэтому центральная двойная связь видна (по межатомным расстояниям) довольно четко.
Сейчас попробую свести в кучку заряды и порядки связей (в т.ч. из своих старых расчетов аналогичных соединений)

[sanya1024]

Базис 6-31G(d,p). Я ничего не пересчитывала, только подняла свои старые выдачи, поэтому для релаксированного S0 цифры из PBE0, а для релаксированного S1 -- из B3LYP (если не ошибаюсь, lmihael пользуется B3LYP -- или нет?). Заряды левдиновские -- в данном случае по ним лучше видно.

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

Заряд на N
Релакс. S0 (PBE0)         Релакс. S1 (B3LYP)        
S0                        S0                        
аминогруппа -0.066308     аминогруппа -0.048076   
гетероцикл   0.032731     гетероцикл   0.044926   
S1                        S1                        
аминогруппа  0.135435     аминогруппа  0.167236   
гетероцикл  -0.050352     гетероцикл  -0.061287   

Порядки связей
Релакс. S0 (PBE0)         Релакс. S1 (B3LYP)        
S0                        S0                        
C-C(Ph)  1.251            C-C(Ph)  1.206   
C=C      1.508            C=C      1.535   
C-C(Het) 1.266            C-C(Het) 1.233
S1                        S1                        
C-C(Ph)  1.207            C-C(Ph)  1.146   
C=C      1.475            C=C      1.532   
C-C(Het) 1.248            C-C(Het) 1.205

Таким образом, по зарядам и порядкам связей нет особой разницы между точками, полученными оптимизацией геометрии S0 и S1. Что вполне логично: судя по форме эксп. полосы поглощения стириловых красителей, сдвиг минимума в S1 относительно S0 невелик.
В состоянии S0 заряд (формально +1, сидящий на азоте гетероцикла) более-менее равномерно размазался по всей молекуле. В состоянии S1 он более-менее собрался на аминогруппе.
Что касается порядков связей, то они в значительной степени задаются геометрией. Релаксированная геометрия S0 -- "хиноид" с тенденцией к выравниванию связей, тогда как в релакс. S1 связи альтернируют. Тут уже некорректно сравнивать данные, полученные с разными функционалами (сходство в зарядах -- оно такое, что его не замазала даже разница в функционалах). Пересчитывать, честно говоря, лень.
Главное, что я хотела сказать -- что полярность состояний при 0 град не меняется при "передвигании костяшек счетов" (изменении геометии, соответствующем релаксированному S0 или S1): атом туда, атом сюда, все сдвиги -- практически вдоль длинной оси молекулы, и по-любому при возбуждении положительный заряд устремляется на аминогруппу.

[sanya1024]

Сделаем корректное сравнение: CIS -- он и в Африке CIS.

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

Заряд на N
Релакс. S0                Релакс. S1        
S0                        S0                       
аминогруппа -0.156380     аминогруппа -0.142810   
гетероцикл  -0.035437     гетероцикл  -0.024691   
S1                        S1                       
аминогруппа -0.097926     аминогруппа -0.079863   
гетероцикл  -0.054519     гетероцикл  -0.054432   

Порядки связей
Релакс. S0                Релакс. S1        
S0                        S0                       
C-C(Ph)  1.217            C-C(Ph)  1.180   
C=C      1.576            C=C      1.626   
C-C(Het) 1.236            C-C(Het) 1.193
S1                        S1                       
C-C(Ph)  1.265            C-C(Ph)  1.234   
C=C      1.387            C=C      1.416   
C-C(Het) 1.266            C-C(Het) 1.237

Снова видно, что в S1 состоянии положительный заряд устремляется с азота гетероцикла на аминогруппу. И снова видно, что порядки связей больше зависят от заданной геометрии (к-ра, к слову сказать, не является в CIS оптимальной ни для S0, ни для S1 -- делались только single-point расчеты). И порядки связей при возбуждении меняются как раз в сторону "хиноидности", т.е., уменьшения разницы между формально двойными и формально одинарными связями. Связано ли это с тем, что с точки зрения CIS мы не находимся ни в минимуме S0, ни в минимуме S1 -- не знаю, проверять лень.

[lmihael]

Всем доброй ночи!

А почему, собственно, эта связь (а также и другие связи на изначальном рисунке) должны так уж легко вращаться?
Если глянуть на дипольные моменты S0 и S1, при возбуждении происходит полная перезарядка, дипольный момент меняет направление. Как мне кажется, у S1 ведущая конфигурация та,
которая отвечет образованию хиноидных структур с заменой одинарных связей на двойные и наоборот. Это - в газе. Что в растворе, сказать намного сложнее.

Здравствуйте Alex. Спасибо что присоединились к обсуждению.

Базис 6-31G(d,p). Я ничего не пересчитывала, только подняла свои старые выдачи, поэтому для релаксированного S0 цифры из PBE0, а для релаксированного S1 -- из B3LYP (если не ошибаюсь, lmihael пользуется B3LYP -- или нет?).

Да, я использовал B3LYP.

Что интересно в древней полуэмпирической статье Cao и др. вращение в газовой фазе было не выгодно по 2-ой связи (вроде ~12 кДж/моль), а в растворителе даже слабополярном (eps=5) картина менялась на противоположную см. рисунок 9. Но у меня такое мощное влияние слабополярного растворителя что-то вызывает сомнения.

Сао.pdf

Может быть мы действительно ищем несуществующую черную кошку в темной комнате? Тогда получается что TDDFT в таких расчетах все таки нельзя доверять?

[sanya1024]

Возбужденное состояние -- с сильным переносом заряда, а следовательно, жутко сольватохромное. На него любой растворитель будет влиять. Но еще и в модели сольватации дело, конечно. Сао и компания использовали модель типа примитивной Онзагеровской, но с цилиндрической, а не сферической полостью, и трехслойной (!) диэл. проницаемостью. От модели сольватации результат может очень даже сильно зависеть (у меня получалось, что разные модели для стириловых красителей давали противоположные тенденции в серии).
Посмотрим на эксп. данные у Сао: сольватохромизм поглощения присутствует (не такой зверский как у Нильского Красного, но заметный), а кв. выход флуоресценции падает в полярной среде на порядок, как и кв. выход транс-цис-фотоизомеризации. Это все -- свидетельства в пользу перехода "хиноид -> бензоид" при возбуждении.
Сольватохромизма флуоресценции особо нет (3-4 нм -- это, считай, ничего), но здорово уширяется полоса поглощения. Большой (>100 нм) стоксов сдвиг. Значит, положение минимума S1 все-таки сильно смещается как по ядерным координатам, так и по энергии.
По их расчетам получается, что для вращения вокруг формально двойной связи при 90 град есть резкий минимум -- могла бы поверить, но только при этом не должно быть щели между S0 и S1. Для остальных вращений в газ. фазе при 90 град есть слабо выраженный минимум. Причем это все -- по картинкам для нерелаксированных сканов. Релаксированные сканы они делали для S0, а для возбужденных считали вертикальные спектры в этих точках. В результате даже те слабо выраженные минимумы просто исчезли.
Single-reference методам в таких расчетах можно доверять только после серьезной проверки multi-reference методами. Но беда в том, что сами multi-reference методы здесь содержат достаточно много произвола: в CASSCF мы можем получить только те возб. состояния, к-рые возможны в рамках нашего акт. пространства. Выкинули одни орбитали, добавили другие -- качественная картинка может измениться. Включили в усреднение не то состояние -- опять качественная картинка изменилась. А QDPT -- совершенно не black-box метод и сильно зависит от качества исходного CASSCF.

[lmihael]

Но беда в том, что сами multi-reference методы здесь содержат достаточно много произвола: в CASSCF мы можем получить только те возб. состояния, к-рые возможны в рамках нашего акт. пространства. Выкинули одни орбитали, добавили другие -- качественная картинка может измениться. Включили в усреднение не то состояние -- опять качественная картинка изменилась. А QDPT -- совершенно не black-box метод и сильно зависит от качества исходного CASSCF.

Но сколько мы не бьемся:
И активное пространство от 2х2 до 16х12 варьировали.
Разные стартовые орбитали из разных точек из разных методов пробовали.
Переставляли орбитали.
А результат практически не меняется - deltaE находится в пределах 14-20 кДж/моль для второго угла.
Может так оно и есть?
Может TDDFT нельзя доверять? А минимум на 90 градусов (если он есть) обусловлен влиянием растворителя?
Тогда остается как-то учесть растворитель на результате QDPT. Весь вопрос как это сделать и как это сделать так чтобы хотя бы самого себя убедить в правильности такого подхода .