Электронный спектр нейтрального радикала

[AllAdd-In]

Пожалуйста, подскажите программу/метод, с помощью которых можно посчитать спектр поглощения хромофора в белке в основном состоянии И в виде нейтрального радикала, то есть окисленном на один электрон (мультиплетность 2) и без одного протона. Хромофор представляет собой ковалентно сшитые боковые группы тирозина и гистидина (к несчастью, гистидин в нем координирован медью). Нужны энергии электронных переходов, их силы осцилляторов и направления дипольных моментов переходов. При этом, разумеется, очень желательно учесть эффекты растворителя (лучше EFP).
Из того, что знаю, методом TDDFT могут считать Firefly, Gamess-Us и Gaussian, но:
1) FF не считает TDDFT для открытых оболочек
2) Gamess-us, насколько я понял, не может учесть растворитель для UHF-TDDFT, а опция spin-flip (для которой это возможно) работает для мультиплетности референсного состояния 3 и больше.
3) Gaussian платный и дорогой.
Насколько я могу судить, остается КВ и CASSCF/XMCQDPT. Если кто считал что-то подобное, поделитесь опытом, плиз.
И кажется, sanya1024 грозилась запостить что-то вроде мануала по расчету спектров. Ох как бы я был благодарен!
Спасибо.

[sanya1024]

Для мануала мне нужен месяц спокойной жизни
Вам подойдет только CASSCF/XMCQDPT в FireFly. Растворитель там, в принципе, можно учесть в PCM или EFP, но в state-averaged варианте, т.е., с одинаковой поляризацией окружения для всех состояний. Если все состояния одинаковой полярности -- это ничего, но если разные (а белковые хромофоры стопудово имеют состояния разной полярности) -- ошибка будет большой. См. тут: http://classic.chem.msu.su/cgi-bin/ceil ... 346-00.htm и дальше. Можно пытаться делать state-specific расчеты (светят проблемы типа вариационного коллапса при расчете возб. состояния) либо state-averaged с игрой в веса, но все это надо тестировать.

Да, CASSCF/XMCQDPT с EFP у меня работало, но у меня хромофор был -- пирен, у него все состояния неполярные.

[AllAdd-In]

Спасибо!
Как только будет время соорудить пособие по спектрам - я первый на очереди!
Если попытаться сформулировать Ваш ответ в привычных мне словах, правильно ли я понимаю, что на данный момент в методе CASSCF/XMCQDPT в FireFly можно учесть растворитель, но его взаимодействие с хромофором будет высчитываться на основании электронной плотности, усредненной по нескольким состояниям (основное и несколько возбужденных)? При этом стационарные дипольные моменты разных состояний могут быть ориентированы по-разному, поэтому усреднение по состояниям может неверно описывать взаимодействие с окружением. Соответственно, у меня есть 2 варианта:
1) попытаться работать без усреднения по состояниям - тогда есть риск, что расчет возбужденных состояний не сойдется
2) попробовать разные веса в усреднении, а правильность описания каким-то образом оценивать по результату. (Интересно, по какому критерию?)
Верно я Вас понял?

[sanya1024]

Да, все совершенно верно. А с п. 2 -- критерий только один, совпадение с экспериментом. Если для Вашей системы эксперимента нет, то для начала потренироваться на похожей, для к-рой эксп. есть. Ну, и не бояться потратить силы и получить отрицательный результат.
Веса можно попробовать взять такие: 1 для интересующего состояния, 0.1 для всех остальных. А программа уже сама нормирует это все к 1. Можно посмотреть, насколько сильно будет меняться результат в зависимости от весов (понятно, что смотреть надо на веса, уже нормированные к 1). Кстати, всем (в т.ч. и Грановскому) будет интересно посмотреть.

[AllAdd-In]

Ну что же, если это интересно не только мне, тогда и вправду стоит попробовать.
Правда, не могу похвастаться богатым успешным опытом расчетов, поэтому буду советоваться и задавать, возможно, глупые вопросы.
Примерный workflow для моей системы ("программа максимум") мне видится следующим образом:
1. составить систему "хромофор, координированный медью 2+, которая, в свою очередь, еще координирована 2мя гистидинами и водой"
2. окружить ее эффективным фрагментом (EF) на свой выбор - из pdb-структуры белка
3. выбрать базис (хм, это само по себе отдельная задача)
4. посчитать single-point на уровне UHF с мультиплетностью 2 и зарядом +2 (все из-за меди)
5. оптимизировать, взяв орбитали из $VEC предыдущего расчета (EF заморозить)
6. посчитать CIS, чтобы определить примерно активное пространство
7. сделать расчет MCCSF с усреднением/без усреднения, причем, насколько я понимаю, CISTEP=ALDET
8. сделать расчет XMCQDPT2, используя оптимизированные орбитали, полученные на пред. стадии
9. поставить мультиплетность 3 (спин 1/2 от меди(II) плюс 1/2 от собственно радикала) и повторить п.п.6-9 (только сначала поставить GUESS=HUCKEL или сделать single-point в оптимизированной геометрии, но для радикала)
как-то так?

экспериментальные спектры поглощения есть только для модельных соединений Гис-Тир (только ковалентно сшитые боковые группы, без protein backbone), без меди, в водном буфере.
наверное, можно в качестве "программы минимум" сначала проделать все это просто с модельным соединением и $PCM SOLVENT=H2O $END,
только, соответственно, мультиплетности взять 1, а потом 2 - для радикала.
Буду рад любым конструктивным комментариям.

[sanya1024]

Ну что же, если это интересно не только мне, тогда и вправду стоит попробовать.
Правда, не могу похвастаться богатым успешным опытом расчетов, поэтому буду советоваться и задавать, возможно, глупые вопросы.

И это правильно

Примерный workflow для моей системы ("программа максимум") мне видится следующим образом:

Будем вносить поправки.

1. составить систему "хромофор, координированный медью 2+, которая, в свою очередь, еще координирована 2мя гистидинами и водой"
2. окружить ее эффективным фрагментом (EF) на свой выбор - из pdb-структуры белка

Сложность первая: фрагменты из аминокислот или молекул воды? Любые фрагменты, кроме воды, можно генерить только в виде EFP2. А с EFP2 невозможна оптимизация геометрии. Вода уже существует в виде (параметризованного) EFP1, с к-рым возможна и оптимизация, и гессианы, и расчеты возб. состояний.

3. выбрать базис (хм, это само по себе отдельная задача)

Выбор, на самом деле, невелик: что влезет, то и берите. Меньше 6-31G(d,p) брать неприлично, а больше... ну, если влезет -- хорошо, а нет -- что ж поделаешь.

4. посчитать single-point на уровне UHF с мультиплетностью 2 и зарядом +2 (все из-за меди)

ROHF или даже RO-DFT (возможно, пригодятся оба варианта). С UHF Вы получите черт знает что (Different Orbitals for Different Spins), к тому же наборов орбиталей будут два (для альфа и для бета-спинов), и Вы не сможете нормально вставить их в CASSCF.

5. оптимизировать, взяв орбитали из $VEC предыдущего расчета (EF заморозить)

Если фрагменты замораживать, тогда действительно все равно, какие они будут, EFP1 или EFP2.

6. посчитать CIS, чтобы определить примерно активное пространство
7. сделать расчет MCCSF с усреднением/без усреднения, причем, насколько я понимаю, CISTEP=ALDET

ALDET лучше распараллелен (если собираетесь считать на кластере) и эффективнее в случае больших акт. пространств. ALDET неизбежен, если будете усреднять по состояниям разной мультиплетности. Если усредняете по состояниям одной мультиплетности, а акт. пространство небольшое -- может лучше сработать GUGA.

8. сделать расчет XMCQDPT2, используя оптимизированные орбитали, полученные на пред. стадии
9. поставить мультиплетность 3 (спин 1/2 от меди(II) плюс 1/2 от собственно радикала) и повторить п.п.6-9 (только сначала поставить GUESS=HUCKEL или сделать single-point в оптимизированной геометрии, но для радикала)
как-то так?

В качестве GUESS-а можно взять вектора от предыдущего CASSCF (с мульт. 2). Это гораздо лучшее стартовое приближение, чем HUCKEL, и не придется делать все пп., а только пп. 7-9. Да, геометрию все равно придется переоптимизировать, но если геометрия не сильно покорежилась (скажем, не стала циклической из линейной), для нового CASSCF будут лучшим стартом вектора из старого CAS-а, чем из HF.

экспериментальные спектры поглощения есть только для модельных соединений Гис-Тир (только ковалентно сшитые боковые группы, без protein backbone), без меди, в водном буфере.

Начните с этой системы.

наверное, можно в качестве "программы минимум" сначала проделать все это просто с модельным соединением и $PCM SOLVENT=H2O $END,
только, соответственно, мультиплетности взять 1, а потом 2 - для радикала.
Буду рад любым конструктивным комментариям.

Да, сначала модельное соединение в PCM, потом можно взять большой кластер воды (~35-40 молекул) в виде EFP1 и соптимизировать либо в HF, либо в DFT. Можно попробовать совместить (небольшое кол-во EFP-шек, чтобы заткнуть водородные связи, + PCM)

[AllAdd-In]

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

Но модельное соединение (не радикальную форму) можно считать в RHF, там меди-то нет?
или сразу в ROHF?

[sanya1024]

Конечно. Форму с закрытой оболочкой -- в RHF или R-DFT. ROHF-а там никакого не получится, даже если его задать в инпуте.

[AllAdd-In]

Ну я же обещал глупые вопросы
Осмелюсь сразу выложить out из single-point RHF, дабы наметанный взгляд выявил недостатки уже сейчас.
Считаю на СК "Чебышёв"

YH_plane_hf_1point.out.zip

[sanya1024]

Ну, во-первых, на такую маленькую молекулу натравливать аж 64 проца -- это жестоко (по отношению к кластеру). Во-вторых, памяти можно и побольше заказать (60 Мслов -- мало, я не помню сколько памяти на нодах Чебышева, но сотню-то уж можно просить). COORD=ZMT -- уже негуманно по отношению к Вам. Если собираетесь оптимизировать геометрию, лучше COORD=UNIQUE и nzvar=<любое ненулевое значение> в $CONTRL и добавьте группу $ZMAT DLC=.TRUE. AUTO=.TRUE. $END . Это означает задание декартовых координат в инпуте, но оптимизацию в автоматически сгенерированных внутренних делокализованных.
Добавьте в $CONTRL d5=.t. -- использовать сферические гармоники вместо декартовых d функций. Избавите себя от проблем с линейной зависимостью. А проблемы уже появились, обртите внимание на
* * * WARNING * * *
------------------------------------------------------------------------------
THE OVERLAP MATRIX HAS 19 EIGENVALUES BELOW 2.1E-04.
THE SMALLEST OF THESE IS 3.10218E-07.
THIS INDICATES A PARTIAL LINEAR DEPENDENCE IN YOUR ATOMIC BASIS.

В принципе, такое сообщение вылетает часто не по делу, особенно в развесистых базисах. Но Вы расстарались И 2dp, и диффузных функций напихали на все атомы. Я не уверена, что это так уж необходимо, особенно в PCM, где окружение играет роль стенки и не дает электрону убежать слишком далеко. Может быть, имеет смысл оставить диффузные или только на периферийных атомах (водородах), или только на электроотрицательных с помощью опции ELNEG в группе $BASIS (из мануала: ELNEG(1)=7,8,9 will add diffuse functions to N, O, and F, but not to other atom types. Works only when DIFFSP is set to .TRUE.).
Имеет смысл последовать совету программы и поставить INTTYP=HONDO ICUT=11 ITOL=30 в $CONTRL -- это по-любому не повредит. Все манипуляции с базисом и способом расчета интегралов лучше проделать раз и навсегда заранее, до того как начнете ставить серьезные расчеты, поскольку они влияют на величину энергии. Начнете менять по ходу дела -- придется все пересчитывать.

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

[AllAdd-In]

Зато с 64 процами и 60 Мсловами влез аккурат в тестовую очередь, где никого нет, и за минутку посчиталось. В других-то очередях и 200 можно запрашивать.
ок, попробую посчитать с Вашими поправками и отчитаюсь

[sanya1024]

Ну, я скорее про то, что и 32 процов было достаточно.

[AllAdd-In]

Получились значения энергии (а.е.)

-568.4077446390 (1) - первый раз
-568.4006025006 (2) - диффузные функции только на N и O
-568.3966342143 (3) - диффузные функции только на H

Так-то с более полным базисом энергия получается ниже.
Частичная линейная зависимость осталась в обоих случаях (2 и 3), но если раньше плохих значений было 19, то теперь их 1.
Присоединяю out из расчета (2)

YH_plane_hf_1point_new_basis.out.zip

[sanya1024]

Так-то с более полным базисом энергия получается ниже.

Кто б сомневался Но на самом деле нас всегда интересуют разности энергий. Если энергия радикала относительно формы с замкнутой оболочкой мало меняется в зависимости от базиса -- можно взять базис поменьше.
К чему я это: пока мы имеем дело с маленьким хромофором в QM/MM, PCM или EFP, все ОК и большой базис вполне влезает. Но если вдруг нам придется "оживить" часть окружения (т.е. перевести ее в QM часть), задача может раздуться, и придется экономить на всем. Поэтому заранее лучше выяснить, чего нам будет стоить такая экономия.

Частичная линейная зависимость осталась в обоих случаях (2 и 3), но если раньше плохих значений было 19, то теперь их 1.

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

Параллельно поставьте оптимизацию геометрии обеих форм в HF (нам нужно хорошее стартовое приближение для CASSCF). Это тоже быстро проскочит.

Кстати, мне показалось, что WALL CLOCK TIME в случае 64 процов было чуточку побольше, чем на 32? или это считалось на других нодах?

[AllAdd-In]

Нет, узлы те же. Но у нас ведь и базис поменьше.
Из интересу посчитал при прошлых параметрах на 32 процах, получилось в 1,5 раза дольше.
Для случая (2) разность энергий рад.-нерад. на 5.8% меньше относительно случая (1), для случая (3) - на 7.0%
Вроде не так много.
На всякий случай - табличка с полученными энергиями и их разностями.

radical_energy.txt

[sanya1024]

Ну ОК. Значит, базис выбираем по вкусу, а если вдруг придется "оживить" часть окружения или добавить медь -- легко сможем урезать слишком большой (пересчитав все разности в меньшем базисе).
Поскольку в CIS расчете у нас скорее всего вылезут первым делом валентные pi-pi* переходы, а LUMO в большом базисе -- это ридберговы орбитали, то для CASSCF либо придется выбирать виртуальные орбитали руками, "стаскивая" их вниз, к границе HOMO-LUMO, либо сгенерить усредненные по нескольким интересным для нас CIS-овским возб. состояниям натуральные орбитали (последний вариант лучше, но несколько муторнее).

[AllAdd-In]

Сколько же всего интересного мне предстоит узнать!
И конечно же, я наврал. Не 5,8% и 7%, а 0,6% и 0,7%, соответственно. То есть разница вообще в пределах процента.
В принципе, наиболее эффективен оказался базис с elneg(1)=7,8: и ошибка маленькая, и считается существенно быстрее.
Пока займусь оптимизацией геометрии в этом базисе.

[sanya1024]

Сколько же всего интересного мне предстоит узнать!
И конечно же, я наврал. Не 5,8% и 7%, а 0,6% и 0,7%, соответственно. То есть разница вообще в пределах процента.
В принципе, наиболее эффективен оказался базис с elneg(1)=7,8: и ошибка маленькая, и считается существенно быстрее.
Пока займусь оптимизацией геометрии в этом базисе.

Да и я тоже молодец: цифры видела, а проверить не догадалась Ну, раз решили остановиться на базисе с elneg, то так тому и бывать

[AllAdd-In]

Оптимизировалось, хотя и не сразу без мнимых частот. В нерадикальном состоянии пришлось сдвинуть разок вдоль нормальной моды (сокращенный out для последней стадии прилагаю), а для радикала - аж три раза.
Для нерадикальной формы прилагаю CIS-out. Буду рад советам по формированию активного пространства для MCSCF.
Для радикальной формы посчитать не могу: "CIS program ... сurrent implementation allows the use of RHF reference wavefunctions only."

YH_plane_hf_opt_offset_trunc.out.zip

YH_plane_CIS.out.zip

[sanya1024]

Для радикала CIS не надо считать. Соптимизируйте в ROHF, этого достаточно.
Ой, заметила: Вы гессиан в PCM считали. Это некорректно делать, отсюда и частоты странные. Т.е., оптимизацию сразу надо было делать без PCM (думаю, переоптимизировать немного времени займет).
Может быть, где-то в недрах кемпорта лежит ссылка на дискуссию в CCL по этому поводу.