photon190573


Кругом квантЫ!

лирические записки о квантовой химии


[sticky post]ну вот, открываемся :)
photon190573
Это место для телег на тему квантовой химии. Но моя аудитория -- не те продвинутые души, кто навскидку определяет кластерные амплитуды для n-кратных возбуждений и выписывает выражения для любого порядка теории возмущений. Я буду писать для тех, кто хочет продвинуться дальше рутинных квантовохимических расчетов в Гау$$иане. Писать буду по-русски, потому что по-английски и без меня довольно написано. Формул будет минимум, постараюсь обойтись вообще без них, если что -- дам ссылку. Азы объяснять не буду, на то есть учебники. Регулярности в постах не обещаю, все зависит от моего настроения :) Осмысленные комменты принимаю с благодарностью. Бессмысленные -- удаляю, ибо нефиг. Оскорбительные для меня или моей аудитории -- удаляю с проклятьями в адрес оставившего их анонимуса.

Апдейт:
Народ, я совершенно не против анонимных постов, если они по делу. Но мне было бы удобнее отслеживать ход дискуссии не по IP автора поста, а по какому-то более удобному идентификатору :) Если у вас нет аккаунта ЖЖ, входите через другие соцсети. На крайняк -- хотя бы представляйтесь. Заранее благодарна.

Все, что Вы хотели знать об интерпретации выдач TDDFT расчетов, но боялись спросить -- продолжение
photon190573
Что еще полезного можно извлечь из выдачи TDDFT расчета? Например, поляризацию перехода. Ее направление совпадает с направлением дипольного момента перехода (transition dipole). Программы выписывают в явном виде transition dipole для каждого перехода, и это можно изобразить в Кемкрафте стрелкой. Особенность Кемкрафта в том, что для рисования стрелки там надо задавать начало и конец вектора. Удобно, чтобы начало лежало где-то в пределах молекулы (центр молекулы или один из атомов). Конец получится, если к координатам начала прибавить x,y,z компоненты момента перехода. Конечно, для растворов поляризация перехода не имеет особого смысла, все равно там молекулы ориентированы как попало. Но для кристаллов или организованных сред типа жидких кристаллов или упорядоченных агрегатов поляризация перехода -- очень даже полезная инфа.

Часто спрашивают, как определить тип перехода (π-π* или n-π* или что там еще, есть ли перенос заряда, и т.п.). Надо смотреть на вид и локализацию задействованных в переходе орбиталей. Точнее, не во всем переходе, а в его доминирующих конфигурациях.
Вернемся к азам. Вспомним про устройство волновой функции метода КВ (CI) и TDDFT. В методе CIS волновая функция
строится в виде линейной комбинации слейтеровских детерминантов, в к-рых одна из занятых (в основном состоянии) спин-орбиталей заменена на вакантную. Такой детерминант называется конфигурацией, в выдаче указано, какую занятую орбиталь заменили на какую вакантную, и выписаны коэффициенты при конфигурациях. Часто бывает, что в этой линейной комбинации можно выделить один или два доминирующих вклада. Вот их-то мы и высматриваем.
Знак коэффициента большого смысла не имеет, но можно (хотя не особо нужно) смотреть на относительные знаки. Например, если мы нашли две конфигурации с большими вкладами, то знаки их коэффициентов могут быть одинаковыми или разными, и если нашлось одно возбужденное состояние, где эти знаки одинаковые, то где-то поблизости должно быть аналогичное (в смысле, с близкими по модулю коэффициентами), но с разными знаками. Коэффициенты нормированы на 1, т.е. сумма их квадратов равна 1. Соответственно, квадрат коэффициента показывает долю данной конфигурации в общей волновой функции возбужденного состояния, вот и весь их физический смысл.
В TDDFT волновая функция (насколько вообще о ней корректно говорить в случае DFT) устроена чуть сложнее, там участвуют не только EXCITATIONS, как в CIS, но и DEEXCITATIONS. Про последние, не вдаваясь в подробности, надо знать только то, что большой (~0.1) коэффициент при DEEXCITATION косвенно указывает на заметный вклад в данный переход двукратных возбуждений (замен сразу двух занятых спин-орбиталей на вакантные). Вообще говоря, в CIS такого быть не может по построению волновой функции возбужденного состояния, а в TDDFT может, поскольку там при построении возбужденного состояния пляшут не от волновой функции, а от электронной плотности и ее отклика на периодическое возмущение (световой волной). Если коэффициент при девозбуждении ~0.1 или больше -- это плохой признак, это означает, что методом TDDFT для данной системы пользоваться нельзя.
Итак, мы нашли одну или две доминирующих в переходе конфигурации (для этого мы открывали выдачу текстовым редактором). Теперь смотрим, какие именно занятые орбитали заменяются на какие вакантные (а для этого нам нужен визуализатор). Думаю, отличить π орбитали от σ и от n химик сможет. Для не-химиков про это придется делать отдельный пост.
Если начальная и конечная орбитали почти не разнесены в пространстве, такой переход называется локальным. Если занятая орбиталь локализована на одном конце молекулы, а вакантная -- на другом (а то и вовсе на другой молекуле), это означает, что в данном переходе электрон должен перескочить с одного конца молекулы на другой, т.е. должен произойти перенос заряда. Это как раз один из признаков артефактного перехода. Эту "разнолокализованность" задействованных в переходе орбиталей можно оценить количественно с помощью т.н. λ-диагностики. Эту самую диагностику умеет выдавать Гамесс и Гауссиан. Чем меньше величина λ, тем больше степень переноса заряда в процессе перехода, тем менее надежно описывается этот переход с помощью данного функционала.
Если конечная орбиталь расположена где-то на периферии молекулы, причем не столько на атомах, сколько между ними -- это ридберговская орбиталь и переход, соответственно, тоже ридберговский. Не принимайте эти определения буквально, наша классификация вообще довольно условная.

Теперь мы можем сформулировать признаки артефактного перехода. Переход подозреваем на артефактность, если (1) у него почти нулевая интенсивность, (2) подозрительно низкая энергия, (3) по виду задействованных орбиталей или по λ-диагностике понятно, что это перенос заряда на большое расстояние. По отдельности признаки не катят: бывают локальные переходы с нулевой интенсивностью (запрещенные по симметрии, например), низкая энергия -- понятие относительное (для одной системы низко -- это в ИК, для другой -- в видимой области), переходы с переносом заряда бывают довольно интенсивными. А вот все три вместе -- повод задуматься. Если к тому же переход резко едет вверх при увеличении веса ХФ обмена в функционале, то это точно артефакт.
Чего я так распинаюсь по поводу артефактных переходов? Ну, приходилось наблюдать, как не слишком грамотные расчетчики приводят в статье энергию перехода более-менее похожую на экспериментальную. А силу осциллятора предусмотрительно не приводят. И переход характеризуют только как, допустим, HOMO->LUMO, а картинки орбиталей не показывают. И функционал у них -- любимый всеми ламерами B3LYP. Часто простая проверка показывает, что они таки получили артефактное состояние, просто оно случайно попало в более-менее нужный диапазон. А "настоящий" переход в расчет вообще не попал, потому что заказывали 3-4 состояния, и все они оказались артефактами.

Продолжение следует

Все, что Вы хотели знать об интерпретации выдач TDDFT расчетов, но боялись спросить
photon190573
По просьбам дорогих кинослушателей расскажу самую-самую базовую вещь, с к-рой начинается жизнь юного фотохимика-теоретика (паспортный возраст тут не при чем, необходимость заниматься расчетами может настигнуть в любой момент).
Вот вдруг случилось так (ага, должно было случиться, сказал бы Боконон) что Вам пришлось заняться расчетом спектров поглощения или вообще люминесценции, а то и -- о ужас! -- Вас бросили на кинетику и фотохимию. Предположим, Вы даже освоили какую-нибудь программу (будем говорить про Гамесс, FireFly, Гауссиан и Орку), составили инпут и запустили счет. А потом получили выдачу. Боже, как в ней разобраться?
Предположим, Вы открыли Вашу выдачу каким-нибудь визуализатором (я предпочитаю Кемкрафт, но сгодится любой, лишь бы понимал выдачи Вашей программы и умел показывать орбитали). Откройте ее еще и текстовым вьюером или редактором.
Показал ли Вам спектр Кемкрафт или Вы нашли в выдаче энергии переходов и силы осцилляторов. В принципе, это все, что Вам нужно... ну, почти. Считайте, что Вам повезло, если первый переход имеет заметно ненулевую силу осциллятора и этот факт согласуется с экспериментом. Но в жизни бывает по-всякому.
Часто бывает, что запихнув в инпут большую развесистую молекулу с кучей функциональных групп и потратив на расчет туеву хучу времени, получаете на выходе нечто странное: несколько (а то и несколько десятков) переходов с энергией в красном или даже инфракрасном диапазоне и почти нулевой силой осциллятора. А заветный переход в видимой или ближней УФ области, с нормальной силой осциллятора (позволяющей зарегистрировать этот переход в спектре) вылезает хорошо если пятым, а то и тридцатым по счету. Поздравляю, Вы наткнулись на "родовой" недостаток метода TDDFT -- занижение энергий переходов с переносом заряда.
Что это значит? Вообще-то в нормальных молекулах заряды просто так с места на место не скачут, это требует затрат энергии. Но все функционалы в DFT, наследующие LDA, имеют внутренний недостаток: они эти энергозатраты недооценивают, а с некоторыми функционалами вообще получается, что перекинуть заряд с места на место в молекуле ничего не стоит. Именно для борьбы с этой проблемой и вводят гибридные функционалы, в к-рых обменный член содержит примесь точного (Хартри-Фоковского) обмена. В BLYP и PBE добавки ХФ обмена нет -- вот и получаем кучу низколежащих состояний с переносом заряда. Начинаем увеличивать вес ХФ обмена (B3LYP 20%, PBE0 25%) -- эти состояния потихоньку уползают наверх и перестают нам мешать, чем больше вес ХФ обмена, тем дальше они уползают. Собс-но, на этом основан мой любимый тест на "артефактность": берем одну и ту же молекулу и считаем ее спектр с функционалами BLYP (0% ХФ обмена), B3LYP (20%), BHHLYP (50%) и HFLYP (100%). "Нормальные", т.е. локальные состояния, при к-рых электрон далеко не убегает, при увеличении веса ХФ обмена тоже едут вверх, но не так сильно. В принципе, можно выбрать некий компромиссный вариант, при к-ром мешающие переходы уже уехали достаточно высоко, а нужные хоть уехали, но не очень. Еще можно использовать т.н. range-separated функционалы, в к-рых вес ХФ обмена плавно увеличивается с межэлектронным расстоянием (LC-BLYP, wPBE, CAM-B3LYP, wB97 и т.п.). В FireFly их, к сожалению, нет, но есть в Гамессе, Гауссиане и Орке. Такие функционалы меньше портят энергию локальных состояний, но эффективно убирают артефактные состояния с переносом заряда.
Я не хочу сказать, что состояний с переносом заряда вообще не бывает. Бывают, особенно в металлокомплексах. Например, ион [MnO4]-- своей окраской обязан CT переходу. Но там интенсивность-то какая!
Можно ли не обращать внимания на артефактные переходы? Если мы просто хотим построить спектр поглощения -- можно, их все равно видно не будет. Но если нам нужен спектр люминесценции или -- того пуще -- кинетика релаксации возбужденного состояния (т.е., фотофизика или фотохимия) -- лучше от них избавиться. Потому что люминесценция идет из низшего возбужденного состояния, и для того, чтобы правильно ее моделировать, это состояние не должно быть артефактным. А если еще и релаксация возбуждения интересует -- тогда точно между основным и нужным возбужденным состоянием (не обязательно низшим) артефактных состояний быть не должно, все промежуточные состояния должны быть "настоящими".

Продолжение следует

Моя статья на "Биомолекуле"
photon190573
Написала статью на конкурс «био/мол/текст» на "Биомолекуле". Прошу любить, жаловать, комментировать: http://biomolecula.ru/content/2058

Про мертвого лосося
photon190573
Историю про фМРТ мертвого лосося (вроде как даже мороженного) я читала лет этак несколько назад. Интерпретация истории была типично журналистская (ага, не люблю я эту братию): взяли ученые мороженную рыбину, сунули в магниторезонансный томограф, показывали ей картинки, рыбина "реагировала". Вывод: ваше фМРТ -- полное говно, человека от мороженной рыбы отличить не может. Я тогда еще подумала, что наверное в процессе эксперимента рыба начала подтаивать и подтекать, а фМРТ реагирует на ток крови, ну или любого водного раствора, ее заменяющего. Вот на оттаивание полуфабриката девайс и реагировал.
Все оказалось гораздо интереснее (http://brights-russia.org/article/it-died-for-science.html). Оказывается, целью эксперимента было доказать, что "реакция" мертвой (вроде бы свежей, не мороженной) рыбы на стимулы -- это результат не вполне корректной статобработки данных (а сырые данные с томографа невозможно интерпретировать, не обработав как следует). Т.е., если всю статобработку провести как надо (типа, не пожмотиться и выкинуть ненадежные точки, тем самым уменьшив объем выборки), то никакой "реакции" не будет, мертвая рыба будет таки мертвой (договорились -- умерла, значит умерла). Работа была опубликована аж в 2009 или 2010 г. в прикольном Journal of Serendipitous and Unexpected Results и даже получила Шнобелевскую премию. И поспособствовала тому, что гораздо больше фМРТ-шников стали серьезнее относиться к статистической обработке данных, и меньше стало "открытий", к-рые придется закрывать.
К чему я это? последнее время мне несколько раз приходилось разным людям растолковывать, что DFT вообще и TDDFT в частности имеют свои границы применимости.
Например, интуитивно кажется, что теорема Купманса будет выполняться в DFT не хуже, а скорее даже лучше, чем в ХФ (мы же включили корреляцию!), и интерпретировать фотоэлектронные и рентгеноэлектронные спектры в терминах Кон-Шэмовских орбиталей можно так же, как и в терминах орбиталей Хартри-Фоковских. Но пока что известны строгие доказательства аналога теоремы Купманса только для верхней занятой орбитали (HOMO), да и те время от времени подвергают сомнению (сорри, сейчас ссылок не дам, я интересовалась этим года 3-4 назад, может быть, дойдут руки собрать самые свежие статьи по этой теме). Для внутренних орбиталей DFT-шный аналог т. Купманса не доказан, так что интерпретация фотоэлектронных и рентгеноэлектронных спектров по DFT -- на свой страх и риск. Обычно таких интерпретаторов спасает то, что рецензенты тоже не в курсе и полагаются на интуитивное представление :)
Еще одни грабли, на к-рые наступают любители применить готовый метод к чему попало -- расчет высоковозбужденных состояний. Касида (например, тут https://arxiv.org/pdf/1108.0611, тут Annu. Rev. Phys. Chem. 2012. 63:287–323 и тут J. Chem. Phys. 108, 4439 (1998)) показывает, что состояния, лежащие выше потенциала ионизации молекулы (да не настоящего, а того, что дает LDA часть функционала, а это намного ниже), получаются совершенно неправильными из-за неправильного асимптотического поведения LDA части потенциала. Вдобавок (хм, я думала, это очевидно!) среди высоковозбужденных состояний полным-полно тех, что получаются возбуждением двух и более электронов. А эти состояния ни в CIS, ни в TDDFT в привычном нам варианте линейного отклика не попадают чисто по устройству метода -- он только для однократных возбуждений. Поэтому расчет 50 возбужденных состояний "потому что мы это можем" -- дело заведомо бессмысленное. Очень хочется таких товарищей спросить: а какова, по-вашему, природа 48-го возбужденного состояния и какому пику в экспериментальном спектре оно соответствует? и чем вы можете это подтвердить? "Совпадение с экспериментом" при этом часто бывает вполне приличное: спектр высоковозбужденных состояий больших молекул (до начала фотоионизации или фотохимического разрушения) довольно плотный, там всегда найдется что-нибудь плюс-минус похожее на "рассчитанную" энергию.
В общем, если кто-нибудь сделает и опубликует обзор типичных ошибок TDDFT-юзеров, это будет работа, достойная не только миллиарда зимбабвийских долларов :)

Несистематический словарь -- немного о DFT
photon190573
LDA -- Local Density Approximation -- приближение локальной плотности или, что то же самое, приближение однородного электронного газа. Любое DFT начинается с LDA и понемногу улучшается: GGA -- Generalized Gradient Approach -- функционал, включающий помимо собс-но эл. плотности еще и ее градиент. Большинство популярных функционалов -- GGA типа. Meta-GGA -- следующая ступенька "лестницы" -- включает еще и вторую производную (лапласиан) эл. плотности. К meta-GGA функционалам относятся TPSS и Миннесотские функционалы.
Hybrid functional -- гибридный функционал, в к-ром обменно-корреляционная часть выглядит так: E_xc(hybrid) = E_c(pure GGA) + a*E_x(HF) + (1-a)*E_x(pure GGA). Т.е., часть DFT-шного обмена заменена на хартри-фоковский. Самый популярный (и бездумно используемый) функционал B3LYP -- как раз такой, в нем вес ХФ обмена 20%. По стоимости расчетов гибриды сопоставимы с ХФ (чистые GGA дешевле, т.к. не требуют расчета двухэлектронных обменных интегралов). Double-hybrid functional -- дважды гибридный функционал, в нем помимо того, что часть обменного вклада представлена ХФ слагаемым, еще и часть корреляционного вклада заменена на MP2-шный: E_xc(double-hybrid) = b*E_c(pure GGA) + (1-b)*E_c(MP2) + a*E_x(HF) + (1-a)*E_x(pure GGA). По стоимости расчета сопоставим c MP2 (требует соответствующего преобразования интегралов).
Range-separated или long-range corrected fucntional -- гибридный функционал с поправкой на дальнодействие. Вес ХФ обмена в нем, в отличие от обычного гибридного (еще называемого global hybrid), плавно меняется в зависимости от межэлектронного расстояния от 0 (как в LC-BLYP) или от 20% (как в CAM-B3LYP) до 100%. Туда же относятся wPBE и wB97. Стоимость расчета -- как с глобальными гибридами.

Продолжение, конечно, следует...

Спин-орбитальное взаимодействие -- надеюсь, что окончание
photon190573
Итак, после тяжелой и продолжительной... зимы, плавно переходящей в следующую осень... брр... весна и половина лета провалились неизвестно куда. Бывает :(

Короче, сейчас будет обещанное окончание рассказа про спин-орбиту и как интерпретировать выдачу спин-орбитального расчета Гамесса.
Выдача эта организована чудовищно кошмарным образомCollapse )

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

Теперь обещанное примечание про одно- и двухэлектронные вклады в спин-орбиту.Collapse )

А мне понравилось :)
photon190573
Originally posted by superhimik at Так себе стишочек,


но пусть коллеги усмехнутся.



Мело, мело по всей земле

Во все пределы.

Учил кристаллы наш студент,

Учил Валера.


Как летом роем мошкара

Летит на пламя,

Слетались шпоры со стола

К оконной раме.


Доцент чертила на доске,

Но бренным тленом

Была для парня красота

Зон Бриллюэна.


И на засаленный конспект

Ложились тени:

Скрещенье граней и углов,

Осей скрещенье.


И падали карандаши

Со стуком на пол.

Он над контрольной каждый раз

В бессильи плакал.


И всё терялось в снежной мгле

Седой и белой.

И Харгиттаи под столом

Лежал без дела.


Он упаковки постигал,

И страх маразма

Его семестр весь накрывал

Разнообразно.


Чертил решётки он Браве,

И то и дело

От ОЦК и ГЦК

Башка гудела.


Методы учета окружения
photon190573
Я помню, что должна была написать, как расшифровывать выдачу спин-орбитального расчета. Но это реально тяжкий труд -- настолько эта выдача запутана. Сейчас же возникла необходимость (собс-но, она давно была, просто руки дошли) написать про различные методы учета окружения. Как-то я писала про довольно экзотический способ работы с окружением, а сейчас -- про наиболее типичные, те, к-рыми пользуюсь сама.

Окружение молекулы можно учитывать либо в модели непрерывной поляризуемой среды (Polarizable Continuum Model, PCM) с главным параметром -- диэлектрической проницаемостью, либо в разных дискретных моделях. PCM хороша тем, что (1) дешевая, (2) поляризуемая, причем самосогласование поляризации среды и исследуемой молекулы стОит недорого и делается на уровне ССП (в MP2 и CC есть свои особенности, я в них не вдавалась, но по крайней мере знаю где посмотреть про это). Недостатки PCM состоят в том, что она (1) параметризуется эмпирически, (2) не учитывает специфическую сольватацию. Т.е., для неполярных растворителей годится, для полярных апротонных -- худо-бедно, по всякому бывает, для протонных фигово. Можно сделать гибрид: насовать несколько живых молекул растворителя и окружить это все поляризуемой средой. Но может выйти тяжеловато, особенно если 1-2 молекулами растворителя не обойтись, да и знать надо, куда их совать -- такие системы не всегда имеют единственный глубокий минимум. А в белках континуальные модели вообще не работают.

Выход был найден в виде QM/MM (Quantum Mechanics/Molecular Mechanics) -- наиболее существенная часть системы (активный центр, хромофор и т.п.) описывается квантовомеханически на любом желаемом уровне, а окружение -- с помощью мол. механики. Плюсы: (1) окружение более-менее дискретное, в нем можно пытаться описывать специфические эффекты, (2) модель остается относительно дешевой. Недостатки: (1) силовые поля параметризуются эмпирически, и результат зависит от использованного силового поля, (2) большинство используемых в QM/MM силовых полей не поляризуемые, а имеют фиксированные заряды, приписанные к ММ атомам, т.е., невозможно самосогласование взаимной поляризации акт. центра и окружения. Более того, в большинстве программ (в частности, в Гамессе) реализован только т.н. механический эмбеддинг (mechanical embedding), т.е. влияние окружения на акт. центр учитывается чисто механически (шарики-пружинки), и даже те несчастные точечные заряды, приписанные силовым полем к ММ атомам, никак не учитываются. Возможно, это и к лучшему: в схеме, где заряды ММ части фиксированы, ошибки от того, что они-то QM часть поляризуют, а сами-то на ее эл. статическое поле никак не реагируют, могут быть больше, чем если вообще выкинуть из схемы электростатику. Недавно вроде бы появились работы по поляризуемым силовым полям (AMOEBA, например), в т.ч. попытки их включения в QM/MM. Практической реализации я пока не видела.

На всякий случай: ONIOM -- это один из вариантов QM/MM, где между QM и MM частью есть еще прокладка в виде полуэмпирического слоя, на мой взгляд, вещь совершенно дурацкая. Сейчас вроде бы чаще делают ONIOM без полуэмпирического слоя, т.е. в чистом виде QM/MM.

И вот EFP -- Effective Fragment Potentials -- это попытка сделать поляризуемое окружение, параметризованное неэмпирически. EFP -- это жесткий фрагмент (целая молекула или ее часть), где все внутренние степени свободы заморожены, поэтому он может двигаться только как целое. Геометрию его можно соптимизировать любым способом, а потом с помощью Гамесса (эта процедура есть только там) сгенерировать EFP. В EFP входят заряды на атомах и некоторых специальных точках (по дефолту это середины связей (bond midpoints), но можно добавить другие), мультипольные моменты вплоть до октупольного, статические поляризуемости, частотно-зависимые поляризуемости, посчитанные для определенного, довольно широкого набора частот, а также специальный локализованный базис и соответствующий фокиан для межмол. взаимодействий. Ну и всякие параметры для демпфирования и экранирования (damping, screening), чтобы задемпфировать на близком расстоянии межмолекулярное взаимодействие фрагментов между собой и фрагментов с QM частью -- потому что на близком расстоянии (~ длины связи) это уже не межмолекулярные взаимодействия, а нечто иное. EFP используются в Гамессе, FireFly и Q-Chem.

Что можно (в принципе) сделать с EFP? можно оптимизировать геометрию молекулы + окружения. Можно даже всю эту конструкцию запихнуть в PCM (внутренняя сольватная оболочка дискретная, все что дальше -- континуальное). Можно гонять мол. динамику или Монте-Карло (глобальную оптимизацию) фрагментов (только для фрагментов, QM части в таких расчетах быть не должно!). Можно считать спектры в поляризуемом окружении (single-point). Можно разрезать сложное окружение, типа белка, на жесткие фрагменты попроще, сделать из них EFP, потом склеить обратно в белок в нужной конформации и использовать эту конструкцию в качестве окружения акт. центра. Но... все это с оговорками.

Сейчас EFP -- штука слегка недоделанная. Возьмем, например, градиенты. Гамильтониан QM-EFP системы устроен так, что в нем есть слагаемые, отвечающие только за QM часть, только за EFP, и за QM-EFP взаимодействие. И градиенты энергии, столь необходимые при оптимизации геометрии, тоже бывают QM, EFP-EFP и QM-EFP. С первыми и вторыми проблем нет, а вот QM-EFP градиенты почему-то вызвали затык.

В чем же дело? EFP в Гамессе бывают двух видов: EFP1 -- полностью параметризованный (по квантовохимическим расчетам), не содержащий ничего, кроме электростатики и всяких параметров межмол. отталкивания, демпфирования и экранирования, и EFP2, генерируемый процедурой MAKEFP -- тот самый, с электростатикой, локализованным базисом, фокианом и прочим. Чтобы получить EFP1, нужно проделать нехилую работу по его параметризации. Ее однажды проделали -- для воды. Для других растворителей -- желающие, welcome :)) Зато такой EFP полностью функционален, для него запрограммированы все градиенты (EFP-EFP и QM-EFP), можно считать спектры, можно оптимизировать возб. состояния -- все супер, но... фактически, только для воды.

С EFP2 все куда сложнее. Хотя его можно сделать для любой молекулы, средств для работы с ним не так уж много. QM-EFP градиенты для EFP2 до сих пор не запрограммированы (хотя уже пару лет как вышла статья с соответствующими формулами). Если в инпут для расчета спектра вставить EFP2 в том виде, как он генерируется Гамессом, то сам же Гамесс и обругается, мол, не умею спектры с EFP2 считать. Проблему, однако, можно обойти, обрезав EFP2 до состояния EFP1: оставить только электростатику (заряды, мультиполи, статические поляризуемости) и скрининг. Все равно при расчете спектра ничего другого от EFP и не требуется, а на абсолютные энергии состояний мы смотреть не будем, только на энергии переходов. Кстати, если у нас есть несколько конфигураций акт. центра с разным расположением QM и EFP частей, то их энергии в таком варианте (с урезанным EFP2) мы тоже сравнивать не можем -- ни в  основном состоянии, ни в возбужденном. Впрочем, в основном состоянии мы можем сделать single-point расчеты с полным EFP2 и сравнить эти энергии.

Градиенты QM-EFP2 авторы метода обещают уже который год :) последнее, что я от них слышала пару месяцев назад -- это "ну вот еще совсем чуть-чуть осталось". Ждем-с :))

Спин-орбитальное взаимодействие -- продолжение
photon190573
Итак, мы поняли, когда нам придется иметь дело с релятивистскими явлениями в квантовой химии. Становится немножко страшно, потому что тяжелые атомы -- они вот они, рядом: Br, I, Sr, Ba, все переходные металлы -- это даже не экзотика типа трансуранов. Неужели для них для всех нужна другая квантовая химия?

К счастью, не всегда. Релятивистская поправка в 1 порядке теории возмущений распадается на 2 слагаемых -- скалярное и спин-орбитальное. Скалярно-релятивистские эффекты (например, известное всем "релятивистское сжатие" в атомах тяжелых элементов) можно учитывать как с помощью специальных гамильтонианов и полноэлектронных базисов (как это сделано в Природе, и есть такая возможность в Гамессе, но не в FireFly), а можно загнать в псевдопотенциал (+ соответствующий базис) и работать дальше как с обычным гамильтонианом. Эта возможность реализована во всех известных мне программах.

Спин-орбитальные эффекты могут проявляться и на легких атомах. Оно смешивает состояния разной мультиплетности и влияет как на энергии состояний, так и на вероятности переходов между ними.

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

Делает он это по теории возмущений. Здесь невозмущенным гамильтонианом H0 служит бесспиновый гамильтониан, а возмущением V - спин-орбитальное взаимодействие (конкретно -- различные варианты спин-орбитального гамильтониана Паули--Брейта). И поскольку мы считаем спин-орбиту между несколькими состояниями, то и гамильтониан нужен включающий несколько состояний, например, гамильтониан CAS-CI. Вычисляются поправки к собственным значениям этого гамильтониана в 1 порядке ТВ и матр. элементы возмущения, т.е., спин-орбитальные каплинги. А параметры расчета задаются в группе $TRANST (и RUNTYP=TRANSITN). Собс-но спин-орбитальный расчет включается ключом OPERAT=HSO-чего-то-там (можно выбрать любой вариант).

Мы должны задать стартовые орбитали (поскольку предполагается CAS-CI, то сами орбитали модифицироваться при расчете не будут, так что лучше чтоб они были "хорошими") в виде группы $VEC1. Почему не просто $VEC, а именно $VEC1? потому что предполагается еще и $VEC2 -- один набор орбиталей для состояний одной мультиплетности, другой -- для другой. Но ограничения, накладываемые на структуру этих VEC-ов таковы, что лучше использовать для обеих мультиплетностей общий набор орбиталей (NUMVEC=1). А именно: неактивные орбитали в обоих наборах должны быть одинаковыми. Но если у нас $VEC1 и $VEC2 происходят из разных расчетов (например, CASSCF отдельно для триплетов и для синглетов), то неактивные орбитали в них, вообще говоря, одинаковыми не будут, и придется руками подставлять блок неактивных орбиталей из одного VEC-а в другой -- манипуляция, чреватая ошибками, к-рые в форматированном файле трудно заметить. Так что проще сделать один общий набор орбиталей (например, сделав CASSCF с усреднением по синглетам и триплетам вместе) и использовать его в качестве $VEC1.

Ключ NUMCI задает количество различных мультиплетностей, для к-рых мы собираемся считать каплинги. Если у нас только синглеты и триплеты, то NUMCI=2. Случаев чтобы NUMCI было больше, мне не попадалось, но не исключаю, что такие бывают. Меньше (NUMCI=1) -- попадалось (тяжелый элемент, большое спин-орбитальное расщепление между состояниями одной мультиплетности).

NFZC -- число занятых неактивных орбиталей, NOCC -- занятые неактивные + все активные. IROOTS -- массив: число включаемых в спин-орбитальный расчет состояний сначала меньшей мультиплетности, потом бОльшей. Т.е., если у нас синглеты и триплеты, то сначала число синглетов, потом число триплетов. Если в нашей системе основное состояние высокоспиновое (триплетный кислород, например, или септет Eu3+), то надо помнить о том, что здесь мы задаем состояния не по возрастанию энергии, а по возрастанию _мультиплетности_. NSTATE -- сколько в каждом CAS-CI расчете состояний диагонализовать. NSTATE по дефолту равно IROOTS, но можно брать с некоторым запасом.

Есть в группе $TRANST и еще один ключик, к-рый требует отдельного объяснения. Это массив ENGYST. Этот параметр (как и многие другие в мануале Гамесса) описан так, что проще посмотреть, как он должен выглядеть, чем пытаться угадать это по мануалу.

В дефолтном варианте (без ENGYST) невозмущенными энергиями (нулевой порядок ТВ) будут собственные значения H0, т.е., энергии состояний CAS-CI. Но мы знаем, что эти энергии недостаточно правильные и, например, в QDPT получаются лучше. Тогда мы включаем ENGYST и в нем задаем энергии состояний из QDPT - сначала для меньшей мультиплетности (синглеты), потом для бОльшей (триплеты). Если синглетов и триплетов не поровну, недостающее добиваем нулями (программа выстроит их в виде матрицы, где по, условно говоря, горизонтали будут синглеты, по вертикали -- триплеты, и эта матрица должна быть квадратной). Теперь программа будет считать невозмущенными значениями числа из ENGYST и вычислять спин-орбитальные поправки к ним. А для каплингов по-прежнему будет использовать собственные состояния H0, т.е., состояния CAS-CI - все равно ничего лучшего у нас нет.

В ENGYST необязательно запихивать значения из QDPT. Возможно, нам больше понравятся цифры из TDDFT, или мы настолько круты, что смогли посчитать энергии в CCSD. Или мы не стали париться и взяли экспериментальные энергии и подставили их. Неважно, программа примет эти значения за чистую монету и выдаст к ним спин-орбитальные поправки. (Поэтому надо быть осторожными с эксп. значениями -- они _уже_ содержат спин-орбиту, а программа-то об этом не знает. Так что в случае ожидаемой большой спин-орбиты в молекуле лучше эксп. цифры не брать).

Независимо от ENGYST, у нас еще есть $VEC1, т.е., орбитали. На этих орбиталях будут построены слетеровские детерминанты, а из них - собственные состояния H0 (см. тут).

Обычно если мы берем ENGYST из QDPT, то у нас уже есть орбитали из CASSCF. Это очень хороший вариант, потому что (1) CAS-CI сойдется моментально, и (2) это будут те же состояния, к к-рым в QDPT вычисляли корреляционные поправки по ТВ во 2-м порядке.

Если мы берем орбитали из ХФ или из DFT, то неизвестно во что их может превратить процедура CAS-CI (это помимо того, что еще и сходиться может долго). Т.е., в ENGYST у нас энергии, к-рые соответствовали состояниям, построенным процедурами CIS, TDDFT или даже CCSD на ХФ или КШ орбиталях, а в спин-орбитальном расчете процедура CAS-CI на тех же орбиталях построит нам вообще говоря _другие_ состояния, и фиг их знает, как они соотнесутся с теми состояниями, энергии к-рых стоят в ENGYST. Если это нас не смущает, то вперед :)

Ну и в группах $DRT1 и $DRT2 мы задаем собс-но параметры CAS-CI (по тому же принципу как мы задавали их для CASSCF) -- $DRT1 для меньшей мультиплетности (синглеты), $DRT2 -- для бОльшей (триплеты).

И еще одна тонкость. Если у нас есть тяжелые атомы с псевдопотенцилом (да, скалярно-релятивистским, конечно), то для компенсации остовных электронов, запихнутых в псевдопотенциал, используют эффективные заряды ZEFF -- свой для каждого элемента. Это некие полуэмпирические параметры, подбираемые под каждый тип псевдопотенциала для каждого элемента (готовые есть только для SBKJC, но в качестве стартовой точки они годятся в т.ч. для моих любимых Штуттгартских ECP). Ну и поскольку уже сделано грубое полуэмпирическое допущение, то не стоит усложнять себе жизнь, и лучше в таком варианте брать OPERAT=HSO1, чтобы не нарушать магический баланс компенсации ошибок :)

Теперь расчет можно ставить. А как интерпретировать выдачу -- расскажу потом.


Спин-орбитальное взаимодействие
photon190573
Давненько я ничего не писала по делу. Ну, так уж случилось, что мне пришлось подробно объяснять пользователю, как считать спин-орбитальные каплинги (spin-orbit coupling, матричные элементы спин-орбитального взаимодействия) в Гамессе. Продублирую и здесь.

Для чего вообще нам нужна спин-орбита? В обычной, нерелятивистской квантовой механике спин вводится как некая дополнительная степень свободы частицы, обладающая свойствами вращательного момента -- в дополнение к обыкновенному механическому моменту количества движения. В нулевом порядке спиновая переменная не связана с остальными координатами электрона -- пространственные волновые функции электрона (орбитали, ага) просто домножаются на спиновую функцию -- получаются спин-орбитали. Но если взглянуть попристальнее, то два вращательных момента -- спин и орбитальный момент -- должны взаимодействовать. Это и есть спин-орбитальное взаимодействие. Релятивистская квантовая механика обосновывает и свойства спина как вращательного момента, и вид спин-орбитального взаимодействия, и многие другие интересные вещи. Необходимую для прочтения методичку А.В. Зайцевского по релятивистской квантовой химии нашла на Кемпорте и тут.
 
Спин-орбитальное взаимодействие обычно мало, но проявляется в особых случаях. Во-первых, релятивистские эффекты проявляются тем сильнее, чем выше скорость движения электрона, величина релятивистских поправок пропорциональна 4-й степени заряда ядра, так что в тяжелых атомах релятивистские эффекты весьма актуальны (да, для 3d элементов тоже). Во-вторых, в окрестности пересечения состояний разной мультиплетности спин-орбитальное взаимодействие будет определять вероятности переходов, запрещенных по спину. Мы помним золотое правило Ферми: в числителе квадрат матр. элемента взаимодействия, в знаменателе -- плотность состояний или (для дискретного спектра) энергетическая щель. Так вот, в окрестности пересечения состояний энергетическая щель будет стремиться к 0, так что каким бы маленьким ни был матр. элемент взаимодействия, за счет знаменателя вероятность запрещенного перехода будет существенно ненулевой.

Так что релятивистские эффекты имеет смысл рассматривать либо когда имеем дело с тяжелыми атомами, либо когда имеем дело с переходами, запрещенными по спину, и окрестностями пересечения ППЭ разной мультиплетности. Чаще всего удается взять два :)

Продолжение следует.

Учебная исследовательская работа
photon190573
Что-то нехорошие ассоциации вызывает у меня название "Учебная исследовательская работа". Вот с этим оно у меня ассоциируется. И еще с этим. С задачами, к-рые уже стопицот раз решены и никакой научной ценности не представляют. С писанными-переписанными под копирку курсовыми и скачанными с инета рефератами. В общем, с любой деятельностью, результаты которой после оценивания, выкидываются на помойку.
Мы в Центре Фотохимии всегда даем нашим студентам реальные задачи. И требуем результатов -- потому что они нам нужны. Некоторые студенты, получив задание собрать литературу по задаче, приходят с пустыми руками: я ничего не нашел. Конечно, не нашел. На сайтах с рефератами про это ничего нет. А у гугля спросить уже слабо. Особенно если спрашивать надо по-английски.
Ну откуда, откуда они такие приходят? Откуда этот школярский подход: дали задачу -- ищи готовое решение? Из школы, вестимо, но неужели за 3 года в институте его не могут вышибить? Нет, как будто еще крепче вколачивают. И ведь эти дети были лучшими учениками в своих школах или хотя бы учились в лучших школах своего города! На каком этапе обучения у них отключается любопытство и включается "сдал--забыл", только еще и мимо мозга?

Кто, интересно, придет к нам в этом году?

Долгое молчание
photon190573
Молчу. Потому что приличные слова кончились давно, а неприличные... неинформативны, пожалуй.

Защита дипломов. 5 курс, специалисты. Студенты, к-рым меняли на ходу учебные планы и перекраивали программу. Путаются в собственных презентациях, плавают на вопросах. Мало того, что они ни черта не знают -- они этого еще и _не_проходили_! Ядреный, мать его, университет. Лучше бы они так теологию учили -- и по ней бы и дипломы писали. Всяко безопаснее для окружающих будет.

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

Через несколько дней будет защита у 4 курса, бакалавров. У этих вообще ни одного нашего спецкурса не было -- а дипломы защищают на нашей кафедре. Видимо, знания у них должны появиться путем божественного откровения.

А преподаватели, вместо того, чтобы учить хоть чему-нибудь этих оболтусов, весь семестр заполняли и переделывали какие-то УМКД и прочие ФОСы. Нигде, нигде и никогда так не мотали мне нервы, как в этом МИФИ на 1/N ставки за 4000 рублей! Пожалуй, уе...ть лопатой надо все Минобразования в целом, а потом каждого чиновника персонально.

Понимаю, что баян, но не могла удержаться :)
photon190573
У нас было 2 версии Гамесса, взломанный ADF, 6 версия Турбомоля, исходники CFOUR и целое множество псевдопотенциалов всех сортов и расцветок, а также Gromacs, Orca, Природа, последний билд Гауссиана и CP2K. Не то, что бы это был необходимый запас для исследования. Но если начал собирать софт, становится трудно остановиться. Единственное что вызывало у меня опасение — это Гауссиан. Ничто в мире не бывает более беспомощным, безответственным и порочным, чем гауссианщики. Я знал, что рано или поздно мы перейдем и на эту дрянь.

Крик души
photon190573
"Поздравляем, теперь ты преподаватель, сотрудник кафедры на 1/N (N стремится к бесконечности) ставки. Ты думаешь, преподавать -- это нести студентам свет знаний? Хрен тебе: твое дело -- писать учебные планы, выдавать какие-то "компетенции", КИМы (что хоть это такое?) и прочую лабуду, чтобы чиновникам от образования было легче тебя оценивать. Учить будешь кого пришлют, и пофиг, что после часа объяснений, чем сигма орбиталь отличается от пи, назавтра они _опять_ этого не помнят. И никуда ты не денешься -- зачет ты им по-любому поставишь. И приготовь ж... под фитили: ты будешь получать их регулярно. Это высшее образование, детка."

Многоконфигурационные методы как дом, который построил Джек
photon190573
В результате общения с пользователем на форуме поняла, что нужно вернуться к азам и объяснить, кто в МКССП на ком стоит.

В начале был атомный базис. Для тех, кто забыл, какого цвета учебник, напоминаю: это набор более-менее водородоподобных орбиталей, центрированных на каждом ядре в молекуле. И пусть базис у нас "в темном чулане хранится / в доме, который построил Джек".

В методе ХФ или DFT волновая функция или, соответственно, эл. плотность описывается с помощью линейных комбинаций этих атомных базисных функций. Для тех, кто все болел, когда это проходили, напоминаю: эти комбинации называются молекулярными орбиталями, а волновая функция (или ее формальный аналог в DFT) представляется в виде слетеровского детерминанта. Коэффициенты в этих линейных комбинациях (т.е., вклад каждой АО) выбирают таким образом, чтобы минимизировать среднее значение гамильтониана с этим слетеровским детерминантом. Полная электронная плотность (а она является наблюдаемой величиной) определяется занятыми орбиталями. Очевидно, что добавление или отрыв электрона существенно меняет распределение эл. плотности в молекуле.

Таким образом, на первой стадии мы имеем слетеровский детерминант, построенный из молекулярных орбиталей, каждая из к-рых представляет собой линейную комбинацию атомных орбиталей ("которые часто воруют пшеницу / которая в темном чулане хранится / в доме, который построил Джек"). Во многих случаях такого описания основного состояния молекулы вполне достаточно.

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

Возбужденное состояние может иметь одну доминирующую конфигурацию, но не обязательно. Так что представление о том, что возбуждение -- это всегда переход электрона с какой-то одной занятой орбитали на какую-то одну вакантную, -- оставьте, пожалуйста, где-нибудь в другом месте.

Этот метод называется методом конфигурационного взаимодействия (КВ, CI) и может быть разного порядка, или кратности, в зависимости от учитываемых возбуждений (CIS, CISD, CISDT и т.д. -- надо расшифровывать или догадаетесь, что речь идет о Singles, Singles+Doubles, Singles+Doubles+Triples и т.д.?) или включать все возможные возбуждения в рамках заданного набора занятых и вакантных орбиталей, называемого активным пространством. Последний вариант называется КВ в полном акт. пространстве (Complete Active Space CI, or CAS-CI). Если акт. пространство включает в себя все занятые и вакантные орбитали, такой вариант CAS-CI называется полным КВ (Full CI).

Итак, теперь у нас есть представление основного и возбужденных состояний в виде линейных комбинаций слетеровских детерминантов разной кратности возбуждения, построенных из молекулярных орбиталей, каждая из к-рых представляет собой линейную комбинацию атомных орбиталей ("которые часто воруют пшеницу / которая в темном чулане хранится/ в доме, который построил Джек").

Метод многоконфигурационного самосогласованного поля (МКССП, MCSCF) в дополнение к оптимизации коэффициентов КВ включает еще и оптимизацию молекулярных орбиталей. Чем больше акт. пространство (т.е., чем сильнее оно приближается к полному пространству МО), тем меньше нужна оптимизация орбиталей. А оптимизация орбиталей есть не что иное как изменение коэффициентов при АО в разложении МО. И это -- чердак в доме старины Джека :) Таким образом, единственная неизменная вещь здесь, фундамент, на к-ром все стоит -- это базис.

Но где же "те два петуха, которые будят того пастуха"? В смысле, где в этой конструкции место усреднению по состояниям?

В принципе, мы могли бы оптимизировать орбитали для каждого состояния в отдельности, и тогда орбитали для каждого состояния получались бы свои. Но единственный критерий при оптимизации МО и коэффициентов КВ -- это минимизация среднего значения гамильтониана, и этот критерий действителен только для основного состояния данной мультиплетности и симметрии. Поэтому, если нам нужно возбужденное состояние той же мультиплетности и симметрии, что и основное, мы не сможем получить его волновую функцию отдельно. Тут-то и придется использовать усреднение по состояниям.

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

Веса состояний никак не оптимизируют. Это не вариационный параметр. Мы их, вообще говоря, от балды задаем. И "эксперты рекомендуют" :) задавать их равными 1 для всех состояний, включенных в усреднение, и 0 для всех остальных состояний. Это будет означать, что нас интересуют все усредняемые состояния в равной мере, и ни одному из них мы не даем предпочтения перед остальными. И это -- единственный более-менее обоснованный подход в выборе весов. И поскольку веса -- это тот параметр, к-рый задаем мы сами, придавать им какой-то физический смысл, интерпретировать -- не надо. Вот не надо, да.

Гамильтониан, чье среднее значение мы минимизируем, сам зависит от матрицы плотности. Поэтому попытки считать CASSCF в разных точках ППЭ с разными наборами и весами усредняемых состояний можно сравнить с попытками посчитать энергию молекулы в одной геометрии с функционалом BLYP, в другой -- с B3LYP, в третьей -- с BHHLYP, и построить по этим данным профиль реакции.

Итак, в сухом остатке у нас следующая конструкция: базис АО --> МО как линейная комбинация АО --> слетеровские детерминанты, построенные из МО --> состояние как линейная комбинация слетеровских детерминантов --> усреднение матрицы плотности по состояниям как способ получить общий набор МО для нескольких состояний и вообще как способ обращаться с возбужденными состояниями. Анализировать вдоль пути реакции можно, например, веса отдельных конфигураций в составе интересующего состояния. Можно (если не в лом) проанализировать даже вклады отдельных АО в МО, к-рая сильнее всех меняется вдоль пути реакции. А вот веса отдельных состояний в усредненной матрице плотности анализировать не имеет смысла, потому что они такие, какими мы их задали, и в процессе расчета их никто не оптимизирует. Выбрать веса из каких-то физических соображений нельзя, т.к. нет таких физических параметров молекулы, к-рые можно было бы сопоставить с весами состояний в ее матрице плотности.

Надеюсь, я смогла изложить достаточно понятно. Спрашивайте, welcome :)

Снова про CASSCF и вообще многоконфигурационные методы
photon190573
Общение с пользователем катализовало меня на пост о целесообразности использования многоконфигурационных методов. В принципе, я об этом уже писала тут, но, видимо, стоит поподробнее. Я просто скопирую из лички свой ответ этому пользователю с некоторой редакцией.

Иногда в литературе попадаются странные вещи. Например, простое одноконфигурационное основное состояние исследуется методом CASSCF или ваще CASPT2. Смотришь на такое чудо мысли и думаешь: автор, зачем ты машину гонял казеную? (ну не на лаптопе же личном, в самом деле, CASPT гонять?) что ты поймать хотел? просто хотел продемонстрировать, что знаешь такое слово -- CASSCF? а поскольку взять акт. пространство приличного размера кишка оказалась тонка (ну или кластер слабоват :) ), то сделали CAS(2,2) -- 2 электрона на 2 орбиталях.
Пользователю, к-рый только начал осваивать многоконфигурационные методы, такое простительно. А автору статьи в журнале просто не попался въедливый рецензент и не заставил объяснять, зачем все-таки многоконфигурационные методы применили к одноконфигурационной задаче.
Итак, когда многоконфигурационные методы нужны, и когда -- нет.

Обычно к многоконфигурационным методам прибегают, когда либо основное состояние подозревается на неодноконфигурационность (открытые оболочки, где N электронов рассаживаются по M квазивырожденным орбиталям я-собьюсь-со-счета-сколькими способами), либо интересует описание основного и пачки возбужденных состояний на одном уровне. Либо исследуют ППЭ, где ожидаются участки с квазивырождением (т.е., фактически сразу придется исследовать целую пачку близколежащих ППЭ). В этих случаях используют усреднение по целевым состояниям.
Использовать многоконфигурационные методы только для основного состояния в случаях без квазивырождения -- стрелять из пушки по воробьям и с большой вероятностью мазать. Если интересует только основное состояние, то усреднение не делают -- а зачем? И тогда получается, что сначала мы в CASSCF слегка улучшим основное состояние добавкой некоторых конфигураций из акт. пространства, а оно обычно достаточно маленькое, и добавка даже десятка конфигураций с весом >=1% (реально там десятка и не наберется, 3-5 штук помимо самой доминирующей, а остальное -- на уровне пыли) погоды не делает, и вообще не факт, что возбуждения именно с этих активных орбиталей существенно улучшат описание корреляции в основном состоянии. И, не удовлетворившись этим, мы дополнительно улучшаем описание основного состояния во 2 порядке ТВ -- делаем расчет QDPT.
В таких случаях показано старое доброе MP2 -- относительно дешево, сердито, и корреляция должным образом учтена, и никакой лишней работы, все в одном флаконе.

Осмысленность совсем уж минималистских акт. пространств (типа (2,2) или (4,4)) тоже, имхо, в основном сомнительная.
Обычно такие акт. пространства берут, когда точно знают, что все квазивырождение обусловлено именно этими электронами именно на этих орбиталях. Например, в случае квазивырождения в открытой оболочке имеет смысл включать в акт. пространство только ее и ничего более. Тогда эл. состояния, возникающие в ней, будут описаны достаточно примитивно -- только в виде необходимых лин. комбинаций конфигураций из этого жалкого акт. пространства, но вся корреляция будет прекрасненько учтена во 2 порядке ТВ.
Другое дело -- когда исследуют возб. состояния. Там обычно априори ничего не ясно. Даже если CIS или TDDFT показывают, что в первом возбужденном состоянии доминирует HOMO->LUMO, какие еще возбуждения дадут вклад в его описание в CASSCF -- хз. Там обычно начинают с маленького акт. пространства (например, (2,2), т.е. HOMO+LUMO, но не обязательно, можно сразу больше брать) и понемножку расширяют. Причем добавлять орбитали по одной (сначала +1 к занятым, потом +1 к виртуальным и т.п.) -- плохая идея, т.к. такие акт. пространства обычно несбалансированы и либо сходимости не будет, либо сойдется к перекошенному решению. Как расширять акт. пространство -- вопрос проб и ошибок и приходит с опытом. Критерием служит заселенность нат. орбиталей (ее показывает Кемкрафт): близкая к 0 или 2 (~0.01 или 1.99, соответственно) заселенность означает, что данная орбиталь вклада в описание усредняемых состояний не дает, и ее спокойно можно выкинуть... иногда. Потому что если исследуется ППЭ, то орбиталь, не важная на одном участке, может оказаться ключевой на другом. Выкинем ее здесь -- не сможем гладко, без разрывов построить кривую "дотуда".

Так что если вдруг в литературе попался CASPT расчет с акт. пространством типа HOMO+LUMO, и при этом ни возб. состояния не интересуют, ни квазивырождения в основном состоянии не предвидится -- перед Вами типичное пускание пыли в глаза. Ну и сами так никогда не делайте :)

CASSCF -- окончание. Анализ выдачи QDPT
photon190573
В прошлый раз я забыла написать про NSTATE в $XMCQDPT. Как я говорила, это размер эффективного гамильтониана.

Что такое эффективные операторы, можно почитать тут: http://www.chem.msu.ru/rus/teaching/zajzevskii/zaizev.pdf (лучше прочитать, чтобы понимать терминологию). Модельное пространство в MCQDPT натянуто на выбранные вектора КВ, полученные (обычно) в результате SA-CASSCF расчета (т.е., как раз те, что нас интересуют: они входили в wstate(1) в CASSCF и входят в wstate(1) и avecoe(1) в QDPT). Эти вектора -- только часть векторов CAS-CI, полное число к-рых равно числу конфигураций в CAS-CI. Если акт пространство небольшое, то полное кол-во конфигураций, соответствующих заданной мультиплетности, будет вполне обозримым. Так, в CAS(2,2) будет всего три синглетных конфигурации (ψA*ψA, ψB*ψB, (ψAα*ψBβ - ψBα*ψAβ)) и одна триплетная (ψA*ψB); в (CAS(4,4) синглетных конфигураций будет 20, а триплетных 15; а в CAS(6,7) септетных конфигураций будет 7. Конечно, в этом случае бессмысленно мудрить и брать размер модельного пространства меньше, чем полное число векторов CAS-CI. Но кол-во состояний в CAS-CI растет очень быстро, угнаться за ним в NSTATE невозможно. Поэтому для больших (>50 конфигураций) акт. пространств надо ограничивать размер эффективного гамильтониана каким-то разумным числом. Поскольку после всех преобразований (к-рые могут занять ооооочень долгое время) длительность дальнейшего QDPT расчета пропорциональна кол-ву состояний в NSTATE, то размер эффективного гамильтониана лимитируется либо таймлимитом Вашего кластера, либо разумным желанием не делать лишнюю работу.

Опыт показывает, что энергия первого возбужденного состояния (S1 в органическом красителе, без вырождения по симметрии) перестает заметно меняться при NSTATE=10. Если нужно больше состояний, то NSTATE надо увеличить пропорционально. Реально больше 40 состояний не берут, 20-30 -- вполне нормально.

Еще (по мотивам вопросов на Кемпорте) одно замечание. Если нам нужны состояния разной мультиплетности, и мы их изначально получили в детерминантном (cistep=aldet) CASSCF расчете с pures=.f. (собс-но, так и надо делать, если нам нужен общий набор орбиталей для состояний разной мультиплетности), в QDPT смешать разные мультиплетности не получится. Надо на одном и том же наборе $VEC делать отдельные QDPT расчеты для каждой мультиплетности. Для этого надо задать соответствующий MULT в $CONTRL и pures=.t. в $DET либо задать MULT в $XMCQDPT.

Итак, XMCQDPT расчет закончился. Ищем в выдаче строчку  *** XMC-QDPT2 ENERGIES ***. В табличке после нее будут вполне понятные вещи: энергии (CAS-CI на векторах из $VEC) исходных состояний и их исправленные во втором порядке ТВ энергии. Также нам могут пригодиться моменты переходов (RADIATIVE TRANSITION MOMENT). Для обычного спектра поглощения нужны будут RADIATIVE TRANSITION MOMENT BETWEEN STATE #   1 AND STATE # (все остальные), но в выдаче также присуствуют моменты переходов между всеми состояниями (уж не знаю, насколько они могут пригодиться). Еще могут понадобиться ZERO-ORDER PROPERTIES FOR STATE # (т.е., натуральные орбитали и распределение заряда в нулевом порядке ТВ для каждого из состояний), в частности, TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS (малликеновские заселенности).

Кое в чем нам поможет ChemCraft. Открыв им выдачу, мы можем увидеть непосредственно электронный спектр, где энергии переходов выдаются в нм, см-1 или эВ, а силы осциллятора вычислены из моментов переходов по известной формуле. Очень удобно. Еще ChemCraft покажет нам, переход с каких орбиталей на какие доминирует в каждом из возбужденных состояний. Tools->Orbitals->Render molecular orbitals -- и мы увидим натуральные (активные) орбитали с заселенностями в каждом состоянии. Заселенности, сильно отличающиеся от 0 или 2, как раз и указывают на орбитали, доминирующие в переходе. А распределение зарядов в каждом состоянии придется вытаскивать руками из TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS, поскольку ChemCraft нам тут не помощник.

ChemCraft показывает форму натуральных орбиталей (Zero order QDPT natural orbitals) для каждого из состояний в NSTATE. То, что для каждого из состояний форма одной и той же орбитали оказывается разной, -- нормально. Это в SA-CASSCF (из-за усреднения) орбитали были общими для всех состояний. Теперь мы каждое состояние уточнили индивидуально -- и вид орбиталей для каждого из них может стать своим собственным.

Напомню еще раз, что State 1 в CASSCF и QDPT -- это основное состояние данной (в MULT в $CONTRL) мультиплетности. Поэтому если нам нужны энергии, скажем, триплетов относительно основного синглетного состояния, то в выдаче для триплетов на энергии переходов, выдаваемые ChemCraft-ом (и отсчитанные от State 1 в данной выдаче, т.е. T1), смотреть не надо, а надо вычислить их ручками -- вычесть из них энергию S0 из выдачи для синглетов.

Некоторые тонкости.

Если посмотреть по строчке EIGENVECTORS OF THE EFFECTIVE HAMILTONIAN, то можно увидеть матрицу эффективного гамильтониана в базисе собственных векторов CAS-CI. Это то, во что превращается CAS-овский гамильтониан (в базисе этих векторов он диагональный) после включения возмущения. А возмущением здесь служат 1- и 2-кратные возбуждения во внешнее пространство. Если вид этого эффективного гамильтониана не сильно отличается от диагонального (на диагонали числа, близкие к 1, вне диагонали -- близкие к 0), то у нас вообще идеальная ситуация -- CAS-овский гамильтониан почти не придется уточнять по ТВ. Такого, конечно, не бывает. Поэтому очень хорошей считается ситуация, когда на диагонали стоят числа ~0.8-0.9, а вне диагонали числа ~0.1 и меньше. Но и так бывает редко. В обычном MCQDPT им. Накано большие внедиагональные элементы приводят к Большой Лаже. В XMCQDPT это не страшно, если перемешивание векторов происходит в пределах массивов состояний, задаваемых wstate(1) и avecoe(1) (т.е., в пределах модельного пространства). Если заметно примешиваются состояния, лежащие за пределами wstate(1) и avecoe(1), это повод получше посмотреть на исходный CASSCF и, скорее всего, что-то в нем серьезно скорректировать и пересчитать.

Если посмотреть по строчке EIGENVALUES OF THE NON-SYMMETRIC EFFECTIVE HAMILTONIAN, то там будет массив (вообще говоря, комплексных) чисел. Если мнимая часть (второй столбец) маленькая (отличная от 0 в 6 или более знаке), то все ОК. Если отличие от 0 появляется в знаке 5-м или меньше -- тоже повод пересмотреть исходный CASSCF.

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

CASSCF -- продолжение. Что делать дальше
photon190573

Итак, CASSCF сошелся. У нас есть чудненькие орбитальки, отлично согласующиеся с нашим представлением о системе. Ну или не очень согласующиеся, но не вызывающие резкого отторжения. Давайте посмотрим на энергии состояний, к-рые мы включали в усреднение. Нравятся? скорее всего, не очень. Иногда -- просто ни к черту. А что б вы хотели? это только начало!

Как я раньше уже говорила, ХФ не включает в себя электронную корреляцию вовсе, а МКССП включает только ее статическую часть. Динамическая корреляция (та, что связана с попытками электронов держаться друг от друга подальше) учитывается либо с помощью дополнительного КВ (уже над полученной многоконфигурационной волновой функцией), либо с помощью теории возмущений. На настоящий момент даже ограниченное (Singles, Singles+Doubles) КВ над CASSCF -- это дорого (хотя в отдельных случаях, на маленьких системах, его можно использовать в качестве бенчмарка), а теория возмущений 2 порядка -- вполне доступна.

Read more...Collapse )

XMCQDPT, реализованная в FireFly, представляет собой улучшение теории Накано, обладающее некоторыми полезными свойствами: (1) явная нетривиальная эффективного гамильтониана от размера модельного пространства в любом порядке ТВ; (2) сходимость энергий и др. свойств по отношению к расширению модельного пространства в каждом фиксированном порядке ТВ; (3) эффективный гамильтониан должен зависеть от всего подпространства, натянутого на на выбранные вектора CASCI, а не от конкретного выбора базиса в этом подпространстве; (4) полученные энергии должны быть однозначными, непрерывными и гладкими функциями от атомных координат и др. внешних параметров, с возможными исключениями в виде точек конических пересечений и др. областей случайного вырождения.


Как перейти от CASSCF к XMCQDPT? Во-первых, нужно оставить весь инпут от CASSCF, но в качестве группы $VEC взять OPTIMIZED ORBITALS от сошедшегося CASSCF расчета. Во-вторых, добавить ключ MPLEVL=2 в группу $CONTRL. В третьих, добавить следующие строчки:

 $XMCQDPT INORB=1 nstate=<кол-во состояний в эффективном гамильтониане>                                                 
  edshft=0.02 wstate(1)=<кол-во единичек, совпадающее с wstate в $DET или $GUGDM2>,-0 avecoe(1)=<кол-во единичек, совпадающее с wstate в $DET или $GUGDM2>,-0 $end
 $mcqfit $end

INORB=1 означает "не пересчитывать CASSCF, взять орбитали из $VEC и все настройки CASSCF из $MCSCF и $DET/$DRT, но проверить орбитали на ортогональность".
edshft -- величина т.н. сдвига энергетических знаменателей. Эта техника применяется для избавления от вторгающихся состояний. Вообще говоря, величина этого сдвига влияет на вычисляемые энергии, поэтому чем она меньше, тем лучше. Прелесть XMCQDPT в том, что для этого метода edshft=0.02 вполне достаточно, и трогать это значение не надо, тогда как в CASPT2, например, применяют значения ~0.3-0.4, и это как бы еще один способ "подрулить" энергии (т.е., еще один источник произвола). Присутствие группы $mcqfit $end включает некие полезные численные алгоритмы, значительно облегчающие жизнь.

Теперь существенные вещи касательно железа: QDPT работает только на одной ноде, но в несколько потоков. Это означает, что пускать его на нескольких нодах кластера бессмысленно, практически все время все ноды кроме одной будут простаивать. Поэтому QDPT мы будем пускать на одной ноде со следующими настройками:

В группу $SYSTEM добавим MWORDS=50 masmem=<намного больше>. Надо так распределить память, чтобы почти вся оперативная память досталась мастер-процессу (masmem), понемногу (50 мегаслов) получили бы остальные процессы, и сколько-то осталось бы системе. Понадобятся также строчки

$system mklnp=<кол-во процессоров на ноде> np=<кол-во процессоров на ноде> $end
$smp smppar=.t. $end

Остальные настройки -- такие же, как были в исходном CASSCF.

Продолжение следует. Анализируем результаты.


CASSCF -- продолжение. Анализ результатов
photon190573

Итак, после первого удачного запуска (т.е., когда программа не выкинула наш расчет с ошибкой сразу же, а честно смолотила все и открутила заданные в maxit 20 итераций) надо посмотреть, что получилось. Открываем выдачу в текстовом редакторе (никогда так не делали?) и ищем слово final. Если в строке FINAL MCSCF ENERGY IS стоят нули, а строчкой выше мы видим слова о том, что расчет не сошелся, и орбитали для рестарта живут в PUNCH файле -- это нормально, в первый раз редко бывает по-другому. Нам надо открыть эту же выдачу визуализатором (я предпочитаю ChemCraft) и смотреть одновременно на орбитали в визуализаторе и на конфигурационный состав состояний в выдаче.
В Кемкрафте идем: Tools->Orbitals->Render molecular orbitals. Там будет два набора орбиталей. Первый -- NATURAL (натуральные орбитали), где для остовных (неактивных) приводятся орбитальные энергии, для активных -- заселенности, а виртуальных (неактивных) не будет вовсе. Второй -- OPTIMIZED (оптимизированные, точнее, пока что недооптимизированные орбитали). Там у неактивных остовных и вакантных приводятся орбитальные энергии, а у активных -- нули. Сейчас будем смотреть на вид натуральных орбиталей и на их заселенности. Поскольку у нас было усреднение эл. плотности по набору состояний, то и заселенности тоже получатся усредненные по этому набору. Чем больше заселенность активной орбитали отличается от 2 или 0, тем активнее она задействована в усредняемых состояниях. Натуральные орбитали с заселенностями ~1.99 и 0.01 можно считать совсем неактивными при данном усреднении. Иногда их безболезненно можно выкинуть из акт. пространства и уменьшить его размер. Но если предполагается сканирование ППЭ, то может так получиться, что орбитали, неактивные в одной геометрии, станут очень даже активными в другой. Так что тут надо думать.

Обычно по виду натуральные и оптимизированные орбитали похожи. Конфигурации, выведенные в выдаче после -MCCI- BASED ON OPTIMIZED ORBITALS, относятся именно к оптимизированным орбиталям. Думаю, подробно объяснять, что написано в выдаче, не стоит, если что-то непонятно -- спрашивайте в комментах. Нам важно посмотреть на сами активные орбитали и усредняемые состояния: соответствует ли акт. пространство и орбитальный состав усредняемых состояний нашим ожиданиям. В любом случае, для дальнейшей работы нам понадобятся OPTIMIZED ORBITALS из *.dat файла. Мы вставим их вместо группы $VEC в инпут.
Если вдруг какие-то орбитали из акт. пространства убежали, впустив вместо себя другие, для начала попробуем побороться с этим тупой перестановкой орбиталей (массив iorder в $guess, не забываем включить перестановку ключом norder). Если какие-то состояния убежали из усреднения, можно попробовать опцию ntrack в $mcscf. Важно то, что на первых стадиях мы корректируем наш расчет, направляя его так, чтобы сходился куда нам нужно. При этом, скорее всего, будет достигнута и оптимальная энергия (вариационный принцип).

В чем тут дело: МКССП -- вариационная процедура, и самая низкая энергия отвечает самому правильному решению. Однако это нелинейная процедура, и минимизация может завести нас в неправильный минимум. Наша задача -- направить минимизацию в правильный минимум. Обычно химическая интуиция подсказывает, какой минимум должен быть более правильным. Но иногда мы ошибаемся, и то, что кажется нам правильным, по энергии упорно оказывается выше того, что получается "само собой". Значит, наши представления об этой системе в чем-то не верны, и надо конкретно разбираться.

Итак, для рестарта вставляем новые OPTIMIZED ORBITALS из *.dat файла вместо старой группы $VEC, вносим в инпут коррективы и стартуем заново. Результат анализируем. Если видим, что все хорошо, и итерации ведут себя прилично, сходятся без скачков, то можно увеличить maxit в $mcscf и довести дело до сходимости. И после каждого запуска анализируем орбитали.

Продолжение следует. Что нам делать с результатами?


?

Log in

No account? Create an account