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

[AllAdd-In]

Да, спасибо, дискуссию нашел, правда, без ссылки на CCL: viewtopic.php?f=71&t=69764
Из нее я вынес для себя:
1) в водных растворах для молекулы с OH-группой (а таковая у меня одна есть - в нерадикальном состоянии) PCM вообще не адекватная модель, поскольку сольватация обусловлена водородными связями, которые, естественно, не описываются в PCM. По-хорошему нужно окружать молекулу EFP из нескольких молекул воды, чтобы имитировать водородные связи, а уже всю эту систему погружать в PCM. Но это тоже на потом.
2) Оптимизировать нужно (если не с EFP) в вакууме, так как посчитать правильно гессиан в PCM все равно не получится - при сдвигах координат нужно бы сдвигать и полости PCM, а этого не предусмотрено.
3) Оптимизацию можно считать адекватной, если молекула выглядит "на взгляд" неплохо: ароматические кольца плоские, связи нормальной длины и т.п. Ну и плюс чтобы мнимых частот не было.

Соответственно, дилетантский вопрос: если я хочу посчитать хромофор из белка (не в водном растворе, а вроде как в белке), нужно ли мне вообще оптимизировать структуру, вынутую из pdb, если она и так выглядит адекватно?

Кстати, нормально ли, что с базисом с d5=.t. и с внутренними координатами в punch-файле все орбитали после 401-й состоят в основном из нулей?

А высказывание:

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

- огорчило (интересно все), но пока возразить, пожалуй, нечего.

И спасибо, что заглядываете в мои out-ы и находите ошибки!

[AllAdd-In]

так. оптимизировал в вакууме нерадикал, все равно есть одна мнимая частота (11.44 I) с интенсивностью 0.00447 - она соответствует вращению метильной группы, как и в прошлый раз.
поставил, кстати, purify=.t. в $statpt - получились нулевые интенсивности для трансляций, и ненулевые (порядка 0.01-0.04) для вращений. интересно, почему?

Вот еще забавное наблюдение: при визуализации вращений и трансляций в wxMacMolPlot согласованно сдвигаются все атомы, кроме водородов. А водороды остаются на месте. Это такая фишка программы или у меня что-то не так?

[AllAdd-In]

а, понял, чтобы элиминировать все вращения и трансляции, нужно еще написать $FORCE PURIFY=.TRUE. $END

[AllAdd-In]

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

YH_opt_CIS.zip

[VTur]

Мнимая частота говорит о том, что у Вас квадратичная силовая постоянная отрицательна, т.е. ППЭ в этом месте выпукла вверх, а не вниз. Структура находится вблизи максимума, а не минимума.

Если интересует только полная энергия, то барьер вращения метила около 10 кДж/моль. Значит провретесь не менее, чем на эту величину.

[AllAdd-In]

хм. Но если она отрицательна, а другие положительны, не значит это, что мы в седловой точке (но не в максимуме)?

[VTur]

В седловой для всей ППЭ, для конкретного сечения вдоль двугранного угла - в максимуме.

[AllAdd-In]

ок, спасибо

[AllAdd-In]

кажется, я эту частоту победил.
вручную развернул метил и пересчитал. Завтра выложу структуру и CIS

[sanya1024]

1) в водных растворах для молекулы с OH-группой (а таковая у меня одна есть - в нерадикальном состоянии) PCM вообще не адекватная модель, поскольку сольватация обусловлена водородными связями, которые, естественно, не описываются в PCM. По-хорошему нужно окружать молекулу EFP из нескольких молекул воды, чтобы имитировать водородные связи, а уже всю эту систему погружать в PCM. Но это тоже на потом.

Да, обязательно. Либо EFP, либо "живая" вода

2) Оптимизировать нужно (если не с EFP) в вакууме, так как посчитать правильно гессиан в PCM все равно не получится - при сдвигах координат нужно бы сдвигать и полости PCM, а этого не предусмотрено.

Да.

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

Ну... зависит от целей. Если только спектр посчитать -- достаточно дефолтного градиента 1e-4 и гессиана без мнимых частот. Если гессиан еще для чего-то нужен (термохимия, поиск переходных состояний, колебательные и вибронные спектры) -- лучше свести геометрию до градиента 1e-5. А "на взгляд" -- это может быть и 1e-3. Если остановиться на таком градиенте, можно застрять на пологом "склоне" ППЭ, ведущем в минимум, отличный от полученного "на глазок".

Соответственно, дилетантский вопрос: если я хочу посчитать хромофор из белка (не в водном растворе, а вроде как в белке), нужно ли мне вообще оптимизировать структуру, вынутую из pdb, если она и так выглядит адекватно?

Все-таки нужно. Во-первых, в рентгеновской структуре -- посаженные неким стандартным методом водороды. Во-вторых, см. предыдущий абзац насчет градиента 1e-3 (а в структуре, взятой непосредственно из базы, градиент может быть и ~0.1 -- изменение геометрии на первых шагах будет заметно на глаз)

Кстати, нормально ли, что с базисом с d5=.t. и с внутренними координатами в punch-файле все орбитали после 401-й состоят в основном из нулей?

Да. Это как раз выкинутые орбитали типа XX+YY+ZZ (сферические контаминанты)

вручную развернул метил и пересчитал

Это самый надежный способ

[AllAdd-In]

here it is
(кстати, в параметрах переходов почти ничего не поменялось по сравнению с тем случаем, когда мнимая частота присутствовала)

YH_CIS.zip

[AllAdd-In]

sanya1024, VTur, спасибо за ответы. Многое становится яснее

[sanya1024]

(кстати, в параметрах переходов почти ничего не поменялось по сравнению с тем случаем, когда мнимая частота присутствовала)

А поворот метила на спектре практически не сказывается. Попозже расскажу, как строить CIS-усредненные нат. орбитали для CASSCF.

[AllAdd-In]

Здорово. Жду с нетерпением!

[sanya1024]

Итак, мастер-класс по построению CIS-усредненных натуральных орбиталей.
Нам понадобится *.dat (или PUNCH) файл от CIS расчета. Правда, рекомендуется исходный CIS-овский расчет запускать с вот такой строчкой (чтобы фаза орбиталей не менялась при последующем считывании):
$moorth nostf=1 nozero=1 tole=0 tolz=0 syms=1 symden=1 $end
Еще рекомендуется добавлять строчку
$CONTRL WIDE=.t. $END
но я так редко делаю, т.к. панч становится жуткого размера.
Создаем инпут для расчета свойств (runtyp=prop) с той же рекомендуемой строчкой $moorth и ключом WIDE и группами $VEC и $CISVEC из панча в конце инпута.
$contrl cityp=cis runtyp=prop $end
----стандартные строчки инпута----
$guess guess=moread norb=<полное кол-во орбиталей в панче> $end
Очень важно считать из $VEC все орбитали, а не только занятые.
$guess asis=.t. $end
Это нужно, чтобы программа никак не меняла считанные орбитали.
В группу $CIS добавляем следующие ключи: cisprp=.t. istate=-<кол-во усредняемых состояний> RDCISV=.t.
И еще будет группа
$cisprp istate=-<кол-во усредняемых состояний> $end
Отрицательное значение istate как раз и означает усреднение по соответствующему числу состояний.

Запускаем этот инпут (по времени -- считывание и запись в файл, т.е. ок. 1 минуты на 1 ноде). CIS-усредненные нат. орбитали будут в группе $VEC нового панча. Указанием на то, что все прошло как надо, будет строчка перед $VEC:
CIS NATURAL ORBITALS, E(CIS)= -230.3757704834
В принципе, уже все.

Но мы можем захотеть (на самом деле, это обязательно) посмотреть, как выглядят наши натуральные орбитали. Для этого мы создаем еще один инпут, опять-таки с runtyp=prop и группой $VEC из свежеполученного панча. Но cityp=cis надо убрать или заменить на none. И добавить строчку
$guess prtmo=.t. $end
В аутпуте будут напечатаны INITIAL GUESS ORBITALS
Они-то нам и нужны, но визуализаторы (в частности, ChemCraft) эту группу не понимают. Поэтому мы вручную, в текстовом редакторе заменяем

----------------------
INITIAL GUESS ORBITALS
----------------------
на
------------
EIGENVECTORS
------------
И открываем визуализатором. Вуаля!
Обычно в CASSCF идут орбитали, оказавшиеся на границе занятых и виртуальных. На орбитальные энергии их внимания не обращайте: это приписанные им хюккелевские (у нас guess=huckel) энергии, они не имеют смысла.

Остается самый главный вопрос: по скольким состояниям усреднять?

[AllAdd-In]

Во-первых, спасибо (думаю, не только от меня)

Во-вторых, забавно, что wxMacMolPlt рисует орбитали только в том случае, если перед ними стоит строка, содержащая хотя бы слово "FINAL" (начало строки "FINAL ENERGY IS ...."):

FINAL
------------
EIGENVECTORS
------------

В остальном вроде бы все получилось по инструкции (для примера поставил istate=-5)

Относительно istate: меня интересуют переходы с наибольшими силами осциллятора, начиная от примерно 150 нм (сильно меньше и не должны получиться, по идее). При взгляде на CI-SINGLES EXCITATION ENERGIES складывается ощущение, что первых 16ти возбужденных состояний мне хватит с головой. Значит ли это, что и усреднять я должен по 16ти состояниям?

[sanya1024]

Ох ни фига себе! 16 состояний -- имхо, многовато. Неужели действительно нужна область вакуумного УФ (150 нм)? впрочем, 150 нм -- это по-CIS-овски, по-человечески это будет ок. 200 нм, что больше похоже на правду. И судя, например, по выдаче YH_plane_CIS_vacopt2, 150 нм -- это не 16-е, а только 8-е состояние.
Смотря что нужно. Одно дело -- задача о pi->pi* фотовозбужденных состояниях, другое -- о разрыве sigma-связи, например.
Где-то в недрах Кемпорта лежит статья о том, как выбирать орбитали для CASSCF. Мне кажется, тут автоматическая процедура не пройдет, надо ручками.
Конечно, в орбитальном составе страшная каша, доминирующие возбуждения вытащить невозможно, всего по чуть-чуть. В основном состояния до 8-го включительно -- это pi->pi*, а вот 3-е, 9-е -- ридберги в изрядной степени (не зря мы включали диффузные функции). Если интересует фотоионизация, то ридберги нам нужны, процесс идет через них. Но это ведь вакуумный расчет? в окружении возможны 2 варианта: в PCM или EFP окружение просто не даст заряду уйти с молекулы (энергия ридберговских переходов сильно увеличится), а в "живом" окружении заряд может уйти на орбитали молекул окружения (т.н. CTTS -- charge transfer to solvent переходы).
Смотрим выдачу CIS с PCM -- разницы вообще никакой практически. Похоже, что-то не так с этим PCM-ом, потому что там есть переходы с изрядным перемещением эл. плотности (с кольца на кольцо, или с молекулы на фиг), что должно приводить к большому изменению дип. момента и сильной зависимости от растворителя. А у нас вообще ни один из переходов не сдвинулся больше, чем на 1 нм -- непорядок. Надо будет разбираться.
В общем, так: CIS-усредненные берем по 8 или 9 состояниям. А в качестве альтернативы выбираем руками: занятые с 43 по 46, а виртуальные... 47-50, 52, 54, 55, 58. Многовато виртуальных, может быть тяжело для CASSCF. Ну, прогоните 20 CASSCF итераций в таком акт. пространстве, а там видно будет, может и уменьшим.

[AllAdd-In]

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

Смотрим выдачу CIS с PCM -- разницы вообще никакой практически.

Тут какая-то путаница.
У меня же все расчеты в PCMe, кроме оптимизации/гессиана.
Если посчитать CIS в вакууме, то все переходы сдвинуты в длинноволновую область, причем чем выше по энергии, тем больше расхождение (вплоть до 2-3 нм).

Где-то в недрах Кемпорта лежит статья о том, как выбирать орбитали для CASSCF.

Эта? http://onlinelibrary.wiley.com/doi/10.1 ... 8/abstract

спасибо, прочту.
Пока попробую посчитать как Вы посоветовали, а параллельно еще сам тоже подумаю )

[sanya1024]

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

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

Тут какая-то путаница.
У меня же все расчеты в PCMe, кроме оптимизации/гессиана.
Если посчитать CIS в вакууме, то все переходы сдвинуты в длинноволновую область, причем чем выше по энергии, тем больше расхождение (вплоть до 2-3 нм).

Ой. Sorry. Действительно, во всех файлах было PCM -- почему-то я в одном сначала не заметила Но все равно, 2-3 нм здесь -- это вообще не разница.

Эта? http://onlinelibrary.wiley.com/doi/10.1 ... 8/abstract

спасибо, прочту.
Пока попробую посчитать как Вы посоветовали, а параллельно еще сам тоже подумаю )

Она самая.
Но если Вам нужна поляризуемость, то Вам вообще спектр считать не нужно, и CASSCF тоже не нужно.
Может, Вам вот это (из мануала GAMESS-US):

Код: Выделить всё

$CPHF group     (relevant for analytic RUNTYP=HESSIAN)  
 
    This group controls the solution of the response equations, also known as coupled Hartree-Fock. 
 
POLAR = a flag to request computation of the static polarizability, alpha.  Because this property needs 3 additional response vectors, beyond those needed for the hessian, the default is to skip the property.  (default = .FALSE.) 

И еще (оттуда же):

Код: Выделить всё

$LOCAL  

The following keywords relate to the computation of imaginary frequency dynamic polarizabilities.  This is useful in the development of the dispersion energy formula in the EFP2 model, but may also be computed separately, if wished. 
 
POLDYN = a flag to compute imaginary frequency dependent  dynamic polarizabilities (alpha), by analytic means (default=.FALSE., but .TRUE. if RUNTYP=MAKEFP) 

[AllAdd-In]

хм... Я пока не совсем ясно представляю, что буду делать с этими поляризуемостями, если их посчитаю. А со спектром - знаю! Я же направляю осциллирующие диполи вдоль дипольных моментов переходов, а из параметров всех осцилляторов в модели (сила осциллятора, энергия перехода, ширина линии, направление) вычисляю общий тензор поляризуемости (комплексный) и вращательный параметр. К тому же в случае спектров у меня есть хоть какие-то экспериментальные данные.
Но ОК. Допустим, что мы можем полностью пренебречь диссипативной частью для ароматики, поскольку она поглощает при 200-300 нм, а мы смотрим при 400-500 нм. Тогда я мог бы попробовать посчитать зависимость поляризуемости от частоты (длины волны света). Для этого, насколько я представляю, мне может потребоваться RUNTYP=TDHF и $TDHF nfreq=N freq(1)=f1,f2,...fN $END
Но в каком виде программа выдаст мне эти поляризуемости? Будет ли из вывода программы следовать, как направлен осциллирующий диполь? И, насколько я понял, это опять-таки работает только для SCFTYP=RHF (и в FF, и в Gamess-US) , то есть радикал я так не посчитаю.
Честно говоря, я бы пока не забрасывал расчет спектра - в нем, по крайней мере, я знаю примерно, чего хочу получить и каких цифр ожидать.