Спин-орбитальное взаимодействие -- продолжение

[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) -- попадалось (тяжелый элемент, большое спин-орбитальное расщепление между состояниями одной мультиплетности).

*UPD. NUMCI=3 используется, когда есть взаимодействие трех мультиплетностей, в частности, в ионах Ln3+. Например, в ионе Eu3+ основное состояние ⁷F_j, излучающее -- ⁵D₀. Но для аккуратного воспроизведения положения уровней квинтетов (не только ⁵D, но даже и этих) учета спин-орбиты только между квинтетами и септетами недостаточно. Квинтеты (в частности, ⁵D) очень сильно взаимодействуют с триплетами (³P), и это необходимо учитывать в расчете. Только тогда положение излучающего уровня ⁵D₀ воспроизводится правильно.

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, чтобы не нарушать магический баланс компенсации ошибок :)

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