Как назвать тип перехода?

[АлександрI]

Уважаемые форумчане, помогите разобраться с перераспределением электронной плотности при расчёте энергии возбуждённого состояния методом TD-DFT. Переход HOMO 94 - LUMO 95 получился с SAPS 0.97 и -0.06. Прилагаю изоповерхности, полученные по Chemcraft с соответственными названиями. В первом случае электронная плотность сосредоточена на трицикле, а во втором - размазана по всей молекуле, включая заместитель, повёрнутый на 66⁰. В TD-DFT спектре получилась интенсивная палка... Что я могу ещё сказать? (А.С.П.+П.И.Ч.)

[sanya1024]

Если все правильно, то переход (HOMO 94 - LUMO 95 с коэф. 0.97) у Вас в основном локализованный на трицикле, pi-pi* с небольшим весом переноса заряда на повернутый заместитель. Корректнее было бы (правильно тут отдельные специалисты говорят) построить натуральные орбитали для этого перехода, но и то, что Вы построили, вполне сойдет. Если есть другая ситуация (другой расчет или другое возб. состояние), где вес такого HOMO - LUMO перехода маленький, найдите пару орбиталей, для к-рых коэф. больше 0.5, и постройте их. Если получилось, что орбитали более-менее равномерно размазаны по всей молекуле (и такой переход доминирует), то этот переход будет называться делокализованным.
Кстати, в органике "интенсивная палка" характерна скорее для локализованных переходов (но всякое бывает, да)

[alxyppv]

Уважаемые форумчане, помогите разобраться с перераспределением электронной плотности при расчёте энергии возбуждённого состояния методом TD-DFT.

Для визуализации перехода удобно использовать разность электронных плотностей основного и возбужденного состояния. В такой разностной плотности Вы сразу увидите откуда электрон "ушел" и куда "пришел". Т.е. если рассматривать экситон как образование дырки и электрона, то разностная плотность Вам покажет, где локализованы эти самые дырки и электроны. Если Вы работаете с Firefly, то для построения cube файла с разностной плотности надо указать, для какого состояния ее строить (ISTATE=ХХ), ну и обычные ключи для создания cube-файла. Фрагмент входного файла приведен ниже.

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

 $TDDFT NSTATE=24 ISTATE=18 $end
 $ELDENS IEDEN=1 spind=.f. $END
 $CUBE CUBE=.T. nxcube=201 nycube=201 nzcube=201
 x0cube=-16.5 y0cube=-16.5 z0cube=-16.5
 xincub=0.165 yincub=0.165 zincub=0.165 $END

В этом примере проводится рассчет 24 возб. состояний, для 18-го печатается разностная плотность, при этом cube строится с дотсаточно высоким разрешением (201 х 201 х 201). Можно заменить этот $CUBE на более простой вариант:

$CUBE CUBE=.T. mesh=fine $END

но в итоге разрешение будет не очень высоким

[sanya1024]

Алексей, такой вопрос: можно ли для построения cube использовать PUNCH от уже сошедшегося расчета? для $CISVEC этот номер проходит, но насчет $TDVEC я не уверена (хотя никакой принципиальной разницы для программы между $CISVEC и $TDVEC быть не должно).

[alxyppv]

Не совсем уверен, что понял Ваш вопрос. Вы имеете в виду, можно взять вектора и самостоятельно построить разностную плотность? Не знаю, не пробовал - когда мне надо было, я просто "читал" cube из PUNCH. Или Вы имеете в виду, можно ли использовать вектора $TDVEC для нового расчета TD-DFT? По идее - да, но у меня был не совсем удачный опыт (правда очень давно) связанный с наличием симметрии и "вращением" векторов. Но это скорее всего исправлено в новых версиях (или какой-то ключ надо ставить). Не могу сказать что часто этим пользуюсь, так что возможно мои представления об этом уже устарели.

[sanya1024]

Предположим, мы уже сделали расчет TDDFT, но $ELDENS и $CUBE не задали. Теперь у нас есть $TDVEC в панче. Можем ли мы считать для нового расчета $VEC и $TDVEC, чтобы заново не пересчитывать TDDFT, и построить $ELDENS и $CUBE? Для меня-то не проблема пересчитать заново (я работаю на кластере), но тем, кто считает на PC, это может быть актуально. Я, наверное, сама попробую поэкспериментировать -- вдруг получится.

[alxyppv]

По идее - да, можно, для этого они и сделаны.

[АлександрI]

Корректнее было бы (правильно тут отдельные специалисты говорят) построить натуральные орбитали для этого перехода

Натуральные орбитали есть в панче задачи TDDFT. но как их оттуда извлечь для Chemcraft'a? Это те, о которых идёт речь? Если да, то как их оттуда извлечь для Chemcraft'a? Ведь впереди есть другие группы $VEC, которые и обрабатывает визуализатор...

[sanya1024]

Точно, натуральные орбитали как раз нашла -- внимательно изучила панч. Вырезать их ручками (ой!) или скриптик написать, да чтоб нужный кусок выкусывал... засунуть в инпут как $VEC, guess=moread, а дальше... надо уточнить, на сайте у Грановского есть рецепты для чего-то похожего. Короче, смысл в том, что эти орбитали надо скушать теперь уже из инпута и перепечатать как есть (ничего больше не делая) в аутпут. И тогда уже ChemCraft-ом смотреть.

[АлександрI]

$TDDFT NSTATE=24 ISTATE=18 $end
$ELDENS IEDEN=1 spind=.f. $END
$CUBE CUBE=.T. nxcube=201 nycube=201 nzcube=201
x0cube=-16.5 y0cube=-16.5 z0cube=-16.5
xincub=0.165 yincub=0.165 zincub=0.165 $END
В этом примере проводится рассчет 24 возб. состояний, для 18-го печатается разностная плотность,

Спасибо за ответ. Если я правильно понимаю этот фрагмент, при, например, NSTATE=4 и NSTAT=1 будет печататься (я ранее CUBE не применял, объясните - в каком месте панча - в $CUBE?, которые надо далее вырезать для Chemcraft'а) разностная плотность для 1-го возб. состояния? И вообще, сделав рестарты с NSTAT=M будут получаться панчи с данными для построением Chemcraft'ом М изоповерхностей разностной плотности?

[alxyppv]

Если я правильно понимаю этот фрагмент, при, например, NSTATE=4 и NSTAT=1 будет печататься (я ранее CUBE не применял, объясните - в каком месте панча - в $CUBE?, которые надо далее вырезать для Chemcraft'а) разностная плотность для 1-го возб. состояния?

да, правильно понимаете, NSTAT=1 (это, кстати, значение по умолчанию).
Из панча надо "вырезать" кусок, начинающийcя с $CUBE и заканчивающийся ближайшим $END. В сформированной cube файле строчки с $CUBE и $END надо удалить.

И вообще, сделав рестарты с NSTAT=M будут получаться панчи с данными для построением Chemcraft'ом М изоповерхностей разностной плотности?

да, правильно, но возможно даже рестарты делать не надо. Достаточно давно я обсуждал этот вопрос с А.Г., и он сказал, что сделать так, чтобы в одном расчее печатались все нужные переходы, не очень сложно. Вопрос в том, дошли ли у него до этого руки. Если да, то даже и рестарты делать не надо, ну а если нет - тогда как Вы сказали, делать рестарты меняя номер требуемого состояния. Надо бы этот вопрос на форуме Firefly задать...

[АлександрI]

натуральные орбитали как раз нашла -- внимательно изучила панч. Вырезать их ручками (ой!)

Вырезал без ой...

засунуть в инпут как $VEC, guess=moread, а дальше...

Засунул...
А пока жду результаты, попробовал нарисовать молекулярную диаграмму с числами рядом с атомами, равными разности т.н. зарядов из аутпута для основного и первого возб. состояния. Тогда, в соответствии с указанием Alxyppv, рассматриваю экситон как образование дырки и электрона - и видно: откуда электрон ушёл (отрицательные числа), куда пришёл (положительные числа). Результат: 1:1 соответствует (качественно) изоповерхностям по ChemCraft'у!
Ложка (большая) дёгтя: о зарядах из аутпута FireFly приято говорить не иначе, как с негодованием... Да, это лишь способ охарактеризовать электронную плотность в окрестностях атома сферой общепринятого радиуса. Но опытные люди мне подсказали, что использование зарядов по Бейдеру более обосновано... Подскажите, как их посчитать (без G, у меня FF)

[alxyppv]

Ложка (большая) дёгтя: о зарядах из аутпута FireFly приято говорить не иначе, как с негодованием... Да, это лишь способ охарактеризовать электронную плотность в окрестностях атома сферой общепринятого радиуса. Но опытные люди мне подсказали, что использование зарядов по Бейдеру более обосновано...

Даже не знаю как этот фрагмент комментировать. Проблема не в FF, а в методе расчета зарядов по Малликену (и их Вам выдаст любой код). Это заряды могут получиться какими угодно, особенно если Вы работатет с хорошими базисами (когда базисы небольшие, то заряды получаются разумными, но мне очень не нравится такой подход, когда для получения нужного результата мы играем с базисом). В то же время, я не уверен что разность Малликеновских зарядов для основного и возб. состояной будет так уж плоха. Это надо проверять. И откуда Вы взяли этот бред про сферу общепринятого радиуса?

Но опытные люди мне подсказали, что использование зарядов по Бейдеру более обосновано... Подскажите, как их посчитать (без G, у меня FF)

Как бы я это делал:

1) посчитать плотность основоного состояния и сделать cube файл с достаточно хорошей сеткой (не менее 200 х 200 х 200)
2) загрузить программу Bader (http://theory.cm.utexas.edu/bader/, бесплатная, работатет с кубами). С ее помощью определить границы Бэйдеровских атомов. Заряды, которые получатся при численном интегрировании будут скорее всего плохими из-за остова, но Вас абсолютные заряды здесь мало интересуют (если интересуют - посчитайте AIMAll)
3) сделайте cube с разностной плотностью инресующего Вас состояния с той же геометрией и такой же сеткой, что и файл плотности основного состояния
4) проинтегрируйте его с учетом границ атомов, которые Вы определили в пункте 2. Поскольку остова в разностной плотности нет, численное интегрирование по сетке даст Вам вполне приемлимый результат (это легко оценить - сумма разностных зарядов должна дать ноль). ВАЖНО: мы предполагаем что границы атомов при возбуждении меняются слабо.


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

Если предположение о слабом изменении границ атомов при возбуждении Вас не устраивает, можно сделать по другому: сложить разностную плотность и потность основного состояния. В итоге получится плотность возбужденного состояния. Далее можно посчитать заряды Bader' для полной плотности . Но при этом будут проблемы с интегрированием остова. азность зарядов можент все равно получитья разумной, ибо ошибки сократятся (опять же, можно проверить просуммировав разностные заряды - должен получиться ноль)

как-то так

[АлександрI]

Большое спасибо за ответ.

я не уверен что разность Малликеновских зарядов для основного и возб. состояной будет так уж плоха.

Ура! Моя молекулярная диаграмма имеет право на жизнь

И откуда Вы взяли этот бред про сферу общепринятого радиуса?

Не сердитесь, это я сам такое надумал...
Сейчас я почти весь в борьбе с бейдеровской программой
Рассматривая визуализацию орбиталей в основном состоянии, я заметил, что на изоповерхности можно различить на трицикле два идентичных хромофора из пирокатехиновых ядер, практически не сопряжённых ни друг с другом, ни с п-нитрофенилом (что естественно), через "центральные" атомы кислорода и кислорода (АО из аутпута векторов):

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

119  O   7  Z  -0.001112   
123  O   7  Z  -0.001653   
127  O   7  Z  -0.001979   
131  O   7  Z  -0.000243   
136  O   7 XZ   0.000063  
137  O   7 YZ  -0.012203 


180  C  10  Z   0.000461 
184  C  10  Z   0.000773 
188  C  10  Z   0.000797 
193  C  10 XZ   0.000059 
194  C  10 YZ  -0.016072 

...разве что чуть-чуть с помощью d-орбиталей.
А в возбуждённом состоянии эл. пл. делокализуется по всей молекуле и получается высокая палка в TD-DFT-спектре для этого перехода.
Вопрос: а может эта палка потому высокая, что на самом деле состоит из двух сливающихся палок, обусловленных каждым из двух хромофоров?
(Заранее прошу прощения за использование современных и классических терминов)...

[АлександрI]

Точно, натуральные орбитали как раз нашла -- внимательно изучила панч. Вырезать их ручками (ой!) или скриптик написать, да чтоб нужный кусок выкусывал... засунуть в инпут как $VEC, guess=moread,

...
Сделал аналогичное этому, и что из аутпута, что из панча получил два одинаковых набора орбиталей. Оба они внешне похожи на исходные, но явлются их зеркальным отражением относительно плоскости трицикла. Это что, так и должно быть?
И как их теперь называть - изоповерхности визуализации натуральных орбиталей?
Далее сравнил обе структуры по Сhemcraft'у и получил RMSD=0.075. Но ведь делалось,

ничего больше не делая

откуда такое отличие? Тем более подозрительно, что я уже получал это число при сравнении структур в вакууме и после D-PCM ... ИМХО: на самом деле это погрешность программы FireFly при расчёте моей молекулы и, как первые, так и обсуждаемые структуры идентичны, что кажется вполне естественным
Сначала я грешил на Сhemcraft, проверил, это не так...
Прав ли я?

[sanya1024]

Мне кажется, что не очень. В аутпуте после TDDFT -- канонические ХФ (КШ) орбитали (я не уверена, что можно заказать печать натуральных орбиталей в аутпуте, надо посмотреть). В панче есть и те, и другие, важно не ошибиться и вырезать правильный кусок. В принципе, Ваш переход достаточно однозначно описывается указанной Вами парой молекулярных орбиталей (о чем говорит коэф. ~0.9). Так что если Вы действительно вставили в инпут натуральные орбитали для MOREAD, то разницы между MO и соответствующими NO может и не быть (с точностью до фазы).

"Ничего не делать" со входными орбиталями для FF означает добавить в инпут соответствующие указания ( $moorth nostf=1 nozero=1 tole=0 tolz=0 syms=1 symden=1 $end чтобы не ортогонализовать входные орбитали перед расчетом и $guess prtmo=.t. asis=.t. $end чтобы печатать входные орбитали, причем без изменения). Потом в аутпут после строчек
----------------------
INITIAL GUESS ORBITALS
----------------------
надо ручками вставить
------------
EIGENVECTORS
------------
чтобы ChemCraft увидел входные вектора. Тогда только Вы сможете увидеть натуральные орбитали.
Так что в данном случае "ничего не делать" означает изрядную ручную работу

Далее сравнил обе структуры по Сhemcraft'у и получил RMSD=0.075.

Тут, пожалуйста, конкретнее, что Вы делали? Обе структуры -- это какую с какой? Сравнивали геометрии? (именно это делает ChemCraft, причем если вдруг поменялся порядок атомов, то полная фигня получается). 0.075 -- многовато для программного нуля. Для геометрий (если все корректно сделано) это различие существенно. Или Вы что-то другое делали?

[АлександрI]

Вы что-то другое делали?

Да, делал. Делал не то, что нужно
Вопрос с цифрой 0.075 снимаю
Спасибо за указания, это для меня очень интересно и важно
Сейчас вот поставил новый расчёт...

[АлександрI]

АлександрI писал(а):
Если я правильно понимаю этот фрагмент, при, например, NSTATE=4 и NSTAT=1 будет печататься (я ранее CUBE не применял, объясните - в каком месте панча - в $CUBE?, которые надо далее вырезать для Chemcraft'а) разностная плотность для 1-го возб. состояния?


да, правильно понимаете, NSTAT=1 (это, кстати, значение по умолчанию).
Из панча надо "вырезать" кусок, начинающийcя с $CUBE и заканчивающийся ближайшим $END. В сформированной cube файле строчки с $CUBE и $END надо удалить.

Сделал это, открыл в СhemCraft'е, получил нечто весьма приятного вида (см. файл вложения). Так это и есть результат визуализации разностной плотности?? И с чем его едят?Объясните, пожалуйста, как ЭТО интерпретировать?
Или я сделал не то, или недоделал?

[sanya1024]

На мой взгляд, это не очень похоже на разностную плотность. Скорее уж на полную. Может быть (я так никогда не делала, поэтому тут только мое предположение), нужно построить такие штуки для основного и возбужденного состояний и вычесть одно из другого. Если я правильно понимаю, строчка $ELDENS IEDEN=1 spind=.f. $END означает расчет полной плотности, а ключик DIFFD=.t. как раз включает разностную (delta density), но разностную кого с кем?

[alxyppv]

Это действительно полная плотность. Я Вам кажется наврал немного... В смысле, скопировал строчку из своего входного файла, но не убедился, что там действительно был расчет разностой плотности. Прошу прощения
Надо вот так
$ELDENS IEDEN=1 spind=.f. diffd=.t. skiphf=.t. $END

разностная плотность задается ключом diffd=.t. Надеюсь, Вы не очень сильно пострадали... А проверять лучше на маленьких системах, а потом уже большие запускать.