photon190573: (Default)
По просьбам населения попробую разложить по полочкам, что такое мультиреференсные методы и как они соотносятся друг с другом.

В HF/DFT волновая функция представлена в виде 1 детерминанта Слейтера. 1 детерминант (или 1 спин-адаптированная комбинация детерминантов) описывает 1 электронную конфигурацию. Соответственно, если состояние (неважно, основное или возбужденное) не описывается 1 конфигурацией (например, частично заполненная оболочка или ситуации, когда состояния примерно одной энергии описываются разными эл. конфигурациями), такое состояние требует описания в виде лин. комбинации спин-адаптированных волновых функций, соответствующих разным конфигурациям.

Например, состояние 2p^1 -- это 3 конфигурации: 2p_x^1, 2p_y^1, 2p_z^1. Если не учесть их все с равными весами, волновая функция состояния будет перекошенной (например, будет только p_x в сферически симметричной задаче, что сходу неверно).

Соответственно, любое описание состояния через лин. комбинацию конфигураций, является многоконфигурационным (да, в CIS или TDDFT возб. состояния тоже являются многоконфигурационными, даже если доминирует одна). И не так уж важно, получено ли это описание вариационным методом (любые методы ССП) или по теории возмущений, или через кластерное разложение.

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

Но Full CI, если это не очень маленькая система -- это непозволительно дорого. Количество конфигураций растет как факториал (быстрее, чем экспонента). Хочется как-то облегчить себе жизнь. А давайте ограничимся только 1-2-кратными возбуждениями по отношению к референсной ("возбуждения" при этом могут иметь ту же энергию что и референс) -- это CISD (КВ1+2). Про CIS тут даже и говорить как-то неприлично... Другой вариант: а давайте учтем все возможные конфигурации, но на ограниченном наборе орбиталей (например, в пределах 1 оболочки). Это КВ в полном активном пространстве (CAS-CI). Полное -- потому что учтены все возможные для данных активных орбиталей расположения электронов. А активные -- это те самые орбитали, к-рые мы для этих целей выбрали. Ну и можно совсем уж урезать осетра и сделать ограниченное КВ в акт. пространстве. Правда, сейчас так никто не делает за бессмысленностью.

Итак, если речь идет про КВ (CI), то орбитали не трогают, только оптимизируют коэф. КВ, и варианты бывают Full-CI, CAS-CI, CISD и CISD в акт. пространстве.

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

Что значит -- подправлять орбитали? Значит, что получив коэффициенты КВ и построив лин. комбинацию конфигураций, подставим ее в ур. Хартри-Фока и снова диагонализуем методом ССП, но с учетом того, что теперь у нас не 1 детерминант, а целая лин. комбинация. Вот тут и всплывает 1й нюанс: если мы, допустим, для нашей системы 2p^1 будем подправлять орбитали только для 1 состояния (к-рое у нас в ХФ получилось 2p_x^1), то оно и будет у нас получаться дальше. Это называется State-Specific (SS-MCSCF) Как было косое решение в симметричной задаче, так и осталось. Отсюда растут ноги у идеи усреднять по состояниям -- усредняется как энергия, так и матрица плотности, и индивидуальные конфигурации уже не играют такой роли (State-Averaged, SA-MCSCF). В случае вырожденных состояний этот прием отлично работает. Но и в случае квазивырождения, и даже вообще невырождения (когда между усредняемыми состояниями есть заметная энерг. щель) тоже очень неплохо работает: в одном расчете получаются сразу несколько состояний с одинаковой погрешностью. Если в SS-MCSCF каждое состояние получалось в отдельном расчете (только низшее состояние имеет смысл, а остальные хоть и получаются, но смысла особого не имеют), с риском свалиться в вариационный коллапс или с необходимостью вводить симметрийные ограничения, то в SA-MCSCF усредняемые состояния сразу ортогональны друг другу, и никакого вариационного коллапса. И симметрию можно отключить, особенно если и изначально молекула не была симметричной.

Итак, результат вариационного расчета называется референсом. Многоконфигурационный референс получается в результате MCSCF, в отличие от одноконфигурационного (из HF или DFT). MCSCF бывает без усреднения по состояниям (SS-MCSCF) или с усреднением (SA-MCSCF), бывает полным в акт. пространстве (CASSCF) или ограниченным (обозначений не встречала, хотя задать такой расчет в Гамессе или FF легко). Варианта CISD с оптимизацией орбиталей не встречала, хотя технически, наверняка, возможно. Но скорее всего, будет нереально дорого. Похоже что фишка многоконфигурационного референса как раз в том, чтобы получить решение в небольшом акт. пространстве, а потом его улучшать дальше.
Далее, имея на руках многоконфигурационный референс, можно (и нужно!) улучшать его еще дальше за счет взаимодействия состояний из акт. пространства с состояниями из т.н. внешнего пространства -- теми, где участвуют орбитали, не входящие в состав активных. Это все методы с приставкой MR: MRCI (обычно в варианте CISD) над референсом из акт. пространства, MRPT и MRCC -- соответственно, теория возмущений и Coupled Clusters на референсе из акт. пространства. Из них только MRCI -- вариационный.
Грабли там заботливо разложены на каждом шагу, но это уже отдельная история.
Размерность и прожорливость каждого метода -- это, наверное, не ко мне. Лимитировать могут разные вещи. Небольшая (в смысле размера базиса) задачка в огромном акт. пространстве порождает чудовищное кол-во конфигураций, отчего затыкается любой традиционный алгоритм (но специально заточенный под такие вещи DMRG, наверное, нет). Не особо большой SA-CASSCF(8e,8o) для задачи в огромном базисе тоже может стать кластеру поперек горла. Считали нормально CASSCF для N состояний, увеличили до N+1 -- и уперлись в нехватку памяти или безумное время счета.
photon190573: (Default)
И чего это я ушами хлопаю? совсем забыла! Конференция им. Грановского https://www.chemport.ru/webconference2021/ состоялась 26-27 апреля, видео уже давно выложено.

photon190573: (Default)
В результате общения с пользователем на форуме поняла, что нужно вернуться к азам и объяснить, кто в МКССП на ком стоит.

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

В методе ХФ или 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 :)
photon190573: (Default)
Общение с пользователем катализовало меня на пост о целесообразности использования многоконфигурационных методов. В принципе, я об этом уже писала тут, но, видимо, стоит поподробнее. Я просто скопирую из лички свой ответ этому пользователю с некоторой редакцией.

Иногда в литературе попадаются странные вещи. Например, простое одноконфигурационное основное состояние исследуется методом 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, и при этом ни возб. состояния не интересуют, ни квазивырождения в основном состоянии не предвидится -- перед Вами типичное пускание пыли в глаза. Ну и сами так никогда не делайте :)
photon190573: (Default)
В прошлый раз я забыла написать про 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.

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

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

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

Read more... )

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.

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

photon190573: (Default)

Итак, после первого удачного запуска (т.е., когда программа не выкинула наш расчет с ошибкой сразу же, а честно смолотила все и открутила заданные в 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 и довести дело до сходимости. И после каждого запуска анализируем орбитали.

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

photon190573: (Default)

Переходим к практике.

В каких случаях нам может понадобиться МКССП? самый очевидный ответ: если мы нарвались на вырожденное основное состояние. Типа f^6 или f^8 в лантанидах. Но это не единственная ситуация. Чаще бывает, что основное состояние у нас вполне одноконфигурационное, но есть пачка квазивырожденных возбужденных (причем иногда -- разной мультиплетности). Или в равновесной геометрии основное состояние одноконфигурационное, а в переходном состоянии реакции -- существенно неодноконфигурационное, причем с вырождением триплета и синглета (бирадикал). Во всех таких случаях надо делать МКССП с усреднением по состояниям (State-Averaged, или SA).

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

Про _те_самые_ программы я говорить больше не буду, т.к. пользуюсь главным образом программами семейства GAMESS: GAMESS-US и FireFly. Для МКССП идеально подходит FireFly, т.к. на настоящий момент только там работают максимально эффективные и лучше всего распараллеленные алгоритмы CASSCF и XMCQDPT (об этом тоже позже).

По дефолту стартовым приближением любого ССП метода в программах семейства GAMESS является GUESS=HUCKEL. Для обычного ХФ или DFT это вполне нормальное приближение, и в неосложненных случаях меньше чем за сотню итераций мы получим вполне хорошее одноконфигурационное решение. Но если с того же стартового приближения начать МКССП расчет -- проблемы гарантированы. Это так же верно, как то, что "раскаленной докрасна кочергой можно обжечься, если держать ее голой рукой слишком долго; или, что если слишком глубоко разрезать палец ножом, то обычно идет кровь; или ... что если выпить слишком много из пузырька с надписью "Яд", то почти наверняка, рано или поздно, почувствуешь недомогание." Поэтому для начала нам нужно получить приличные стартовые орбитали.

Но если у нас система с вырождением -- мы никак не сможем получить стартовые орбитали в нормальном спин-ограниченном ХФ расчете (ROHF, Restricted Open-shell HF). Значит, нам надо получить орбитали для похожей системы, но без вырождения. Например, для лантанида в конфигурации f^6 или f^8 надо посчитать тот же атом в конфигурации f^7 (полузаполненная f оболочка, вырождения нет, все однозначно), сообразно изменив заряд и мультиплетность. Таким образом мы получим некий плацдарм для размещения лишних электронов или дырок в оболочке. Посмотрим на орбитали в каком-нибудь визуализаторе, чтобы убедиться что получили именно то, что хотели.

Если основное состояние у нас вполне одноконфигурационное, а интересуют квазивырожденные возбужденные, то со стартовыми орбиталями все очевидно: их берут из обычного RHF расчета.

Теперь самое время сказать, что делать МКССП расчеты ограниченной кратности в акт. пространстве имеет смысл только тогда, когда оно ну оооооочень большое (больше 16 электронов на 16 орбиталях). Пока вписываемся в (16,16) -- лучше делать CASSCF, т.е. учитывать _все_ возбуждения. В GAMESS и FireFly это включается опцией FORS=.t. (FORS -- Fully Optimized Reaction Space, синоним CAS). Дефолтное значение этого ключа определяется довольно сложным образом, но если выбираете между КВ методами ALDET и GUGA, то по дефолту FORS=.t., а про другие варианты говорить будем исключительно по просьбе трудящихся (я пока не сталкивалась с необходимостью использовать другие варианты).

Теперь надо определиться с акт. пространством и усреднением по состояниям.
Если имеем дело с вырожденным основным состояниям, то все очевидно (в компенсацию за сложность генерации стартовых орбиталей): активное пространство -- это набор вырожденных орбиталей, а усреднять надо по всей пачке вырожденных состояний. Т.е., в массив wstate(1) нужно поставить столько единичек (усредняем всегда с равными весами! программа потом сама отнормирует все веса на 1), какова кратность вырождения. Так, в уже неоднократно упомянутых системах f^6 и f^8 кратность вырождения 7.
Если нас интересуют квазивырожденные возбужденные состояния, то в компенсацию за простоту генерации стартовых орбиталей нам придется повозиться с определением активных орбиталей и усредняемыми состояниями. Для этого делаем CIS или TDDFT расчет (в зависимости от того, каким методом делали стартовые орбитали, ХФ или DFT, соответственно) и смотрим, сколько возбужденных состояний нам нужно и с каких орбиталей на какие идут эти возбуждения (ес-но, нужны только доминирующие вклады). Возможно, понадобятся отдельные расчеты для триплетов и синглетов. Соответствующие орбитали включаем в акт. пространство (не забываем про ограничение (16,16)!), а в усреднение -- нужную нам пачку возбужденных + основное. Т.е., если делаем расчет только для синглетов, то включаем S0 + пачку Sn, если для триплетов, то пачку Tn, если без разбора спина (и так тоже бывает), то S0 + пачку Tn + пачку Sk (n вообще говоря не равно k).
Для усреднения по состояниям очень важно, чтобы не было перекосов. Так, если есть пачка (квази)вырожденных возбужденных состояний, то в усреднение надо включать ее целиком, ничего не потеряв и не добавив лишнего. Возможно, для этого придется проделать пробный расчет, проанализировать результат и скорректировать wstate. Если вдруг так получилось, что между верхним включенным в усреднение и нижним не включенным оказалась маленькая щель (~0.01 а.е.), это нехорошо. Надо либо верхнее усредняемое выкинуть (посмотрев на щель между ним и нижележащим), либо нижнее не попавшее включить (опять-таки посмотрев на щель).

В группе $MCSCF нам надо задать всякие вещи типа кол-ва итераций и инструмента управления сходимостью, а также задать тип КВ: ALDET (КВ на детерминантах, Ames Lab. DETerminant full CI) или GUGA (КВ на правильных по спину конфигурациях, Graphical Unitary Group Approach). У обоих методов есть свои достоинства и недостатки, можно будет об этом рассказать поподробнее, но не сейчас.

Итак, что нам нужно для CASSCF расчета:
1. Стартовые орбитали ($GUESS guess=moread norb=<полное кол-во орбиталей в $VEC> (возможно, понадобятся norder=1 iorder()) $END
2. Соответственно, группа $VEC, содержащая полный набор орбиталей (при рестарте это будут OPTIMIZED ORBITALS из *.dat файла)
3. Группа $MCSCF, в к-рой нужно указать cistep=aldet или guga, а также задать кол-во итераций МКССП maxit.
4. Группа $DRT для cistep=guga или группа $DET для cistep=aldet. Здесь задается акт. пространство и кол-во итераций КВ itermx. В $DET также задается кол-во состояний nstate, состояние iroot, для к-рого выводить свойства, и массив wstate(1) для усреднения. Для GUGA расчета понадобятся также группы $GUGDIA, $GUGDM и $GUGDM2, где задаются те же nstate, iroot и wstate(1) -- здесь они почему-то оказались в разных группах, так исторически сложилось :)
5. Крайне полезные опции, без к-рых CASSCF расчет будет безнадежно долгим и прожорливым:
  $ciinp castrf=.t. $end
и
  $gugem pack2=.t. $end
В общем-то все, остальное -- тонкости.

По моему опыту, расчеты МКССП не нужно сразу ставить на много итераций (maxit в $MCSCF), лучше сначала несколько раз прогнать с maxit=20, а в особо тяжелых случаях даже с maxit=15 или 10, и после каждого расчета внимательно изучать орбитали, возможно, переставлять что-то с помощью iorder в $GUESS и стартовать с $VEC от последнего МКССП расчета.
Инструменты управления сходимостью FOCAS и SOSCF в $MCSCF действуют альтернативно: или -- или. FOCAS более медленный, но более стабильный, что важно в тяжелых случаях, когда стартовые орбитали очень далеки от оптимальных. SOSCF сходится быстрее, но может пойти вразнос. Есть полезная опция NUMFO, задающая кол-во итераций FOCAS, к-рые нужно проделать прежде, чем перейти к SOSCF.

Продолжение следует. Дальше мы поговорим об анализе результатов CASSCF расчета и что с ними делать дальше.

CASSCF

Nov. 17th, 2013 09:54 pm
photon190573: (Default)
Сколько ж это времени прошло? МАМАДОРОГАЯ! ну ладно, я вернулась. Потому что надоело по 100 раз повторять на форуме одно и то же. Пусть будет глава будущей методички :)

MCSCF (русская аббревиатура МКССП) -- Multi-Configurational Self-Consistent Field -- квантовохимический метод, к-рый используют для генерирования стартовых состояний в случаях, когда обычные одноконфигурационные методы (ХФ, DFT) не годятся.
Хочу подчеркнуть, что MCSCF, как и просто SCF -- это только начало. В обычном SCF электронной корреляции нет вовсе, в MCSCF она т.н. статическая, или нединамическая. По определению, это как раз та корреляция, к-рая существенна для систем, для корректного описания к-рых нужно больше одного (квази)вырожденного детерминанта. Т.е., как раз та, что получается в МКССП. Отличное определение :) собс-но, формальное определение электронной корреляции вообще не лучше.
Давайте попробуем вообразить, как выглядит электронная корреляция. Где-то в этом блоге я уже приводила аналогию с сотрудниками офиса, к-рые терпеть друг друга не могут, но вынуждены работать в одной комнате -- вот и стараются держаться друг от друга подальше. Такое взаимное избегание электронов -- это электронная корреляция вообще. А ее статическая часть -- это когда у сотрудников-электронов есть несколько равноценных по энергии способов рассесться по своим рабочим местам, не мешая друг другу. Проблема выбора -- это страшная штука, и в рамках одноконфигурационного приближения выбрать лучший из равных способов абсолютно невозможно. На это нарывались, к примеру, те, кто пытался сделать ССП расчет для ионов лантанидов с конфигурацией f^6 или f^8. Целых 7 равноценных способов разместить в f оболочке лишнюю дырку или лишний электрон -- процедура ХФ просто теряется от такого изобилия! Выход здесь -- получить орбитали, одинаково хорошо (или одинаково плохо) описывающие весь набор (почти) равноценных (т.е., (квази)вырожденных)) состояний с помощью процедуры усреднения эл. плотности по этому набору.
Таким образом, в методе MCSCF с самого начала строят электронную волновую функцию как линейную комбинацию нескольких слейтеровских детерминантов и часто (хотя и не всегда) используют усреднение электронной плотности по набору состояний, оказавшихся квазивырожденными.

Частным, хотя и очень важным (с практической точки зрения -- самым важным для нас) случаем МКССП является CASSCF -- Complete Active Space SCF, т.е. МКССП в полном активном пространстве. От МКССП вообще оно отличается тем, что выделяют некий набор активных орбиталей (включающий в себя занятые и виртуальные) и учитывают все возможные возбуждения в рамках этого активного пространства. Т.е., со всех занятых на все виртуальные -- полное МКССП в ограниченном акт. пространстве.

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

Уже сама по себе итерационная процедура оптимизации орбиталей -- нелинейная, и в принципе любой ХФ расчет может не сойтись или вообще пойти вразнос, если не включать специальные процедуры (DIIS или SOSCF, демпфирование эл. плотности, сдвиг уровней, смена стартового приближения). Если на это наложить еще и КВ, и сцепить одно с другим -- кто может гарантировать, что расчет вообще когда-нибудь сойдется? Однако тот факт, что молекула существует, означает, что у нее есть некая волновая функция и электронная плотность, а следовательно, существует и некое решение задачи. Т.е., при правильном задании стартовых условий решение найдется. А если упорно не сходится -- надо смотреть, что задано не так.

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

Profile

photon190573: (Default)
photon190573

April 2024

S M T W T F S
 123456
78910111213
14151617181920
2122232425 2627
282930    

Syndicate

RSS Atom

Style Credit

Expand Cut Tags

No cut tags
Page generated Jan. 17th, 2026 08:33 am
Powered by Dreamwidth Studios