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

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

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

Как перейти от 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.

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

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

В каких случаях нам может понадобиться МКССП? самый очевидный ответ: если мы нарвались на вырожденное основное состояние. Типа 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 расчета и что с ними делать дальше.