S1 поверхность (Firefly)

[Reizen]

Для расчетов имеется доступ к университетскому кластеру, но среди ПО нет Firefly (есть Гауссиан03, Gamess, Nwchem и Молпро).

Если вы сейчас говорите о Ростове-на-Дону и ЮФУ - то нет, есть там Firefly.

Другое дело, что "кластер" у нас организован причудливо, и если вам, например, выделят для ваших задач один-два процессора, то это будет ненамного быстрее, чем на персоналке.

[Tonnie Fox]

Для расчетов имеется доступ к университетскому кластеру, но среди ПО нет Firefly (есть Гауссиан03, Gamess, Nwchem и Молпро).

Если вы сейчас говорите о Ростове-на-Дону и ЮФУ - то нет, есть там Firefly.

Другое дело, что "кластер" у нас организован причудливо, и если вам, например, выделят для ваших задач один-два процессора, то это будет ненамного быстрее, чем на персоналке.

С последним Вы в точку попали, насколько я знаю, 2 процессора выделено в целом под хим.фак, при этом расчеты ведут человек 5, и частенько образуется очередь,в которой у меня приоритет равен 0. Так что зачастую быстрее самому провести расчеты.

[VTur]

Не знаю, не знаю. Я считаю там же, но не на химфаке, а на ЮФУ (В.Н. Дацюку мое почтение!). И очередь небольшая и желающих не очень много

[sanya1024]

Ну вот, дошли наконец руки до этой задачи. Раз мы собираемся стартовать с DFT-шных орбиталей, нам понадобятся TDDFT расчеты вертикального спектра в опт. геометрии каждого из участников реакции (реагента, интермедиата и продукта). Посчитайте хотя бы 3 возб. синглета -- нам надо видеть, какие возбуждения в них доминируют, исходя из этого мы будем выбирать акт. пространство. Триплет T1 не всегда имеет такую же электронную конфигурацию, как S1, так что в этом отношении триплеты нам не помощники. Также мне кажется, что действительно, вероятность попасть на триплетную ППЭ при фотовозбуждении намного меньше, чем на синглетную, поэтому те, кто писал, что реакция идет через S1, скорее всего, правы.
Пока, глядя на выдачу с s0, могу предположить, что понадобятся орбитали 20-24 (если тяжело будет, то хотя бы 22-24) и как минимум 25, а возможно еще 26 и 27. Т.е., минимальный вариант -- это орбитали 22-25. Можно начать с этого.

[Tonnie Fox]

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

[sanya1024]

Значит, в бутанале брали (при замкнутой оболочке) 4 занятых и 3 виртуальных МО. Может быть, они были отделены энерг. щелями от соседних и образовывали "естественное" акт. пространство. То, что пару HOMO + LUMO обязательно включать -- очевиднее некуда. Но чем ее дополнить?
Можно по аналогии с бутаналем взять некое "естественное" акт. пространство, т.е. пачку околограничных МО, отделенную щелью от соседей. Но, например, в структуре 1 HOMO (#24) отделена заметной щелью от нижележащих, а между нижележащими каких-то заметных щелей нет. А нам одной 24-й МО явно мало. Так что будем смотреть на конфигурации, доминирующие в возбуждениях. Больше всего нас, конечно, интересует S1, тем более что S2 и S3 лежат намного выше во всех присланных Вами структурах.

Итак, в структуре 1 состояние S1 представляет собой почти чистое HOMO->LUMO возбуждение, а S2 и S3 -- возбуждения с МО #18-23 тоже на LUMO (#25). Но уже в структурах 2 и 3 в S1 состоянии есть вклады возбуждения 24->26. И мне кажется, что не повредит включить 27-ю. Итак, для начала пробуем включить в акт. пространство занятые с #18 по #24 и виртуальные #25-27. Если окажется слишком много конфигураций, крайние МО (18, 19 и 27) выкинем.

Теперь усреднение по состояниям. Это обязательно при оптимизации возб. состояния в CASSCF, иначе при минимизации энергии в ССП (пусть даже МКССП) случится вариационный коллапс -- оптимизируемое возбужденное состояние провалится по энергии ниже основного и задача сломается (не в смысле что программа остановится и обругается, программа может продолжать считать, но выдаваемые цифры будут полным мусором). Усреднение не даст возбужденному состоянию проватиться ниже основного. Для начала попробуем усреднять только по S0 и S1. Обратите внимание: в CASSCF, в отличие от CIS и TDDFT, состояния нумеруются начиная с основного. Т.е., State 1 -- это S0, State 2 -- S1 и т.д. Веса ставим равные, поскольку state-specific оптимизация в state-averaged CASSCF связана с дифференцированием по весу состояния, и что будет при неравных весах -- хз.

CISTEP=GUGA vs. CISTEP=ALDET. На больших акт. пространствах на GUGA может не хватить памяти, тогда ALDET однозначно. Если памяти хватает, то на 1 ноде (я проверяла на 8-ядерной) GUGA чуточку быстрее. ALDET в варианте Грановского лучше распараллелен по нескольким нодам, так что на кластере для больших акт. пространств рекомендуется ставить ALDET с опцией distci=<число ядер>. Но на 1 ноде это неактуально.

Параллельно с CASSCF расчетом в FireFly поставьте оптимизацию геометрии S1 в TDDFT в GAMESS (на кластере). Возможно, полученная структура будет хорошей стартовой геометрией для оптимизации в CASSCF. Обратите внимание: абсолютные энергии DFT, полученные в FireFly и GAMESS, хоть даже с одинаковыми функционалами и базисами, сравнивать нельзя -- для их получения используются по-разному устроенные сетки интегрирования, а абсолютная энергия в DFT зависит от сетки. При этом все разности энергий (энергии переходов, энергии реакций и т.п.) будут одинаковыми (если не забыли в FireFly поставить d5=.t., а в GAMESS ispher=1 -- или наоборот, d5=.f. и ispher=0, соответственно).

[Tonnie Fox]

С помощью примеров с форума Грановского ( http://classic.chem.msu.su/cgi-bin/ceil ... 374+00.htm ), составил инпут для основного состояния системы 1. Все ли корректно?

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

 $CONTRL d5=.t. SCFTYP=MCSCF RUNTYP=optimize MAXIT=100 $END
 $CONTRL NZVAR=3 inttyp=hondo icut=11 $END
 $smp httfix=.f. $end
 $system mwords=50 timlim=600 $end
 $ZMAT DLC=.t. AUTO=.t. $END 
 $basis gbasis=N31 ngauss=6 ndfunc=1 diffsp=.t. ELNEG(1)=7,8,9 $end
 $GUESS GUESS=moread norb=30 $END
 $MCSCF method=dm2 CISTEP=aldet maxit=200 ntrack=2 ISTATE=1 $end
 $DET NCORE=23 NACT=2 NELS=2 NSTATE=4 GROUP=C1 ITERMX=200 PURES=.T. $END
 $STATPT METHOD=GDIIS NSTEP=100 HSSEND=.t. $END
 $DET WSTATE(1)=1,1 cvgtol=1d-7 $END
 $SCF DIRSCF=.t. $END
 $mcscf acurcy=1d-7 ENGTOL=1.0d-12 $end
 $trans cuttrf=1d-13 $end
 $STATPT OPTTOL=1E-4  NSTEP=100  HSSEND=.t. $END
 $FORCE NVIB=2 vibsiz=0.005 $END
 $data
 1_s0
C1
.......................
 $END
 $VEC
.......................
 $END

1. NSTATE=4 поставил потому что со значением 2 обрывались расчеты, видимо, как и в примере Грановского вторым является триплет.
2. Перепробовал кучу комбинаций NACT= NELS=, всегда одна и та же ошибка, кроме случая NACT=2 NELS=2.

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

 **** ERROR, $ALDET INPUT NCORE=  23 NELS=  14 SZ= 0.000
  IS INCONSISTENT WITH $CONTRL INPUT ICH=   0 MULT=   1

3. Для расчета S1 необходимо изменить только значение ISTATE=1 на =2?

[sanya1024]

1. HSSEND=.t. уберите -- слишком тяжело, мало того что каждый геом. шаг -- это 3 CASSCF расчета, так еще таких шагов придется сделать 6N (N -- число атомов). А он к тому же у Вас два раза повторяется.
2. На память не скупитесь, ставьте mwords=260 или сколько Винда позволит, все пригодится.
3. Ошибка -- потому что программа решила, что у Вас NELS=14 при NCORE=23, а в инпуте я вижу NCORE=23, NELS=2 и NACT=2 (решили сделать самый простой CASSCF, (2,2) на HOMO+LUMO). Технически это правильно (хотя я бы все-таки не ограничивалась только этой парой орбиталей). Почему NELS было воспринято неправильно -- хз. Может, стоит попробовать еще раз, тщательно проверив, что все группы начинаются со 2 позиции в строке (строка, начинающаяся с 1 позиции, воспринимается как комментарий и игнорируется).
4. Важные дополнения: $ciinp castrf=.t. $end (если это не так по дефолту)
5. Для начала сведите CASSCF в начальной точке, а maxit лучше берите не слишком большой (<50). Если почему-то ССП не захочет сходиться (запросто!), сможете рестартовать с орбиталей последней итерации, возможно, понадобится их переставить местами. Или захотите добавить в акт. пространство еще пару орбиталей. Лучше поставьте большой ITERMX, чтобы CASCI гарантированно сошлось.
6. Считывать в moread лучше все орбитали. Да, и cvgtol поставьте лучше 1d-8

[Tonnie Fox]

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

          MAXIMUM GRADIENT = 10.8698604    RMS GRADIENT = 4.7084465

      OPTIMIZATION ABORTED.
    -- GRADIENT OUT OF RANGE
    -- MAXIMUM ALLOWED FORCE (FMAXT) =   10.000

пробовал в $STATPT задать DXMAX=0.2 - не помогло

[sanya1024]

Или стартовая геометрия была плохая (с чего стартовали?), или... в общем, лучше выложите полную выдачу с этой ошибкой.

[Tonnie Fox]

спасибо за наводку, воспользовался геометрией, оптимизированной в более широком базисе, и проблема исчезла.

[sanya1024]

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

[Tonnie Fox]

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

Как раз сейчас пытаюсь оптимизировать S0 с активным пространством из орбиталей c 18 по 27, правильно ли задал инпут? В частности интересует, не напортачил ли в $DET

[sanya1024]

Вроде все ок.

[Tonnie Fox]

1. Ну вот, закончил первую серию расчетов для реагента. Экспериментальных данных у меня к сожалению нет по данному соединению, но если сравнивать с бутаналем, то довольно хорошие результаты получились. Если для бутаналя эксп. энергии перехода S0-S1 и S0-T1 составляют 85 и 78 ккал/моль соответственно, то рассчитанные для моей системы - 95,7 и 102,4. То есть энергия S0 занижена на 16 ккал/моль, но разница между S1 и T1 около 7 ккал/моль как в эксперименте. У Дьюара была сходная картина, только точка S0 завышена на 14. Наверное не корректно сравнивать разные системы, но все ж лучше чем ничего.
2. При расчетах S0 всплыли 2 мнимые частоты, как от них избавиться?
3. Сейчас пытаюсь рассчитать триплетный бирадикал, в начале расчета проскакивают строки

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

    Warning - remapping target state #   1 to converged CI state #   2
    Warning - remapping target state #   2 to converged CI state #   1

что они значат?

[sanya1024]

1. Ну вот, закончил первую серию расчетов для реагента. Экспериментальных данных у меня к сожалению нет по данному соединению, но если сравнивать с бутаналем, то довольно хорошие результаты получились. Если для бутаналя эксп. энергии перехода S0-S1 и S0-T1 составляют 85 и 78 ккал/моль соответственно, то рассчитанные для моей системы - 95,7 и 102,4. То есть энергия S0 занижена на 16 ккал/моль, но разница между S1 и T1 около 7 ккал/моль как в эксперименте. У Дьюара была сходная картина, только точка S0 завышена на 14. Наверное не корректно сравнивать разные системы, но все ж лучше чем ничего.

Сравнение не вполне корректное не столько потому, что разные соединения, сколько потому, что наверняка использовалось различное акт. пространство. От этого результат ооочень сильно зависит.
Плюс в небольшом акт. пространстве (а реально нам доступны только такие, и Ваше (14,10) -- тоже небольшое, раз влезло на настольный PC) энергии переходов всегда будут завышены, и скомпенсировать это можно либо КВ над CAS-ом (First-Order CI, Second-Order CI -- одно- и, соответственно, одно- и двукратные возбуждения во внешнее пространство с многоконфигурационного референса), либо -- более экономно -- по теории возмущений (CASPT, MCQDPT и др. разновидности). Но если цель -- исследовать ППЭ возб. состояния и процессы на ней, то важна правильность воспроизведения ее формы (она зависит от типа включенных в акт. пространство орбиталей и лишь во вторую очередь от размера акт. пространства) и ее взаимное расположение относительно ППЭ близлежащих состояний. В данном случае -- с T1, но возможно, что где-то она окажется сближенной и с S2 или с S0.

Смотрю выдачи: на первой же итерации в каждом файле
EXCESSIVE NUMBER OF ITERATIONS...
MCSCF IS NOT CONVERGED!

Никогда не стартуйте сразу оптимизацию. Ваши стартовые орбитали не настолько хороши, чтобы за 40 итераций сразу CASSCF сошелся. Лучше сначала запустите single-point расчет, сведите CASSCF до конца и уже с этих орбиталей стартуйте оптимизацию. В данном случае Вам повезло (хотя на первом шагу CASSCF был далек от сходимости), и на последующих итерациях ничего не пошло вразнос, а спокойно к чему-то сошлось. Но может и не повезти.

Смотрим доминирующие конфигурации. И видим, что и в выдаче для S0, и в выдаче для S1 первое состояние (т.е. по идее S0) имеет доминирующую конфигурацию 2222220200. Вы уверены, что Вам нужно именно это состояние? Как-то это все некрасиво выглядит. Давайте-как возьмем в качестве стартовых геометрий структуры из DFT или ХФ (R-DFT, RO-DFT и TDDFT или RHF, ROHF и CIS, соответственно) расчетов и просто сделаем для них single-point CASSCF в том же акт. пространстве. Хотя бы посмотрим, какие орбитали нам реально нужны, а какие можно безболезненно выкинуть и облегчить себе жизнь.

2. При расчетах S0 всплыли 2 мнимые частоты, как от них избавиться?

Вы все-таки не поленились посчитать гессиан? Смотрим: первая мнимая частота (с бОльшей абс. величиной) -- поворот метильной группы вокруг связи 1-13. Просто руками поверните этот метил на ~60 град. Вторая -- поворот пропильной группы вокруг связи 1-2. Поверните и ее тоже на какой-нибудь заметный угол. Обычно это помогает. Заодно и энергия S0 понизится.

А зачем Вам были нужны эти гессианы, кроме как убедиться, что у Вас минимум (или не минимум)? в T1 и S1, к примеру, я вижу такое же расположение метила и пропила относительно карбонила, как и в S0. Я бы на всякий случай повернула бы их так же, как советую в случае S0. Если оптимизация свалится обратно -- это был минимум. Если найдет другой, поглубже -- значит, так тому и бывать. Но вообще точность расчета численного гессиана невысокая. Так что и доверия этим мнимым частотам или их отсутствию немного. Чтобы численный гессиан был надежнее, надо порог оптимизации снизить хотя бы до 3e-5, а критерии сходимости CASSCF у Вас уже достаточные -- об этом мы заранее позаботились.

3. Сейчас пытаюсь рассчитать триплетный бирадикал, в начале расчета проскакивают строки

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

    Warning - remapping target state #   1 to converged CI state #   2
    Warning - remapping target state #   2 to converged CI state #   1

что они значат?

Этого не пугайтесь, это нормально. У Вас стоит трекинг состояний -- следить за тем, чтобы считать градиент именно того состояния (в смысле конфигурационной природы, а не номера), к-рое Вы выбрали. Просто то состояние, к-рое Вы выбрали в качестве первого, меняется местами со вторым. Если оптимизировали второй триплет, то в результате оптимизации он станет первым.

Но все-таки давайте сначала разберемся с орбиталями. Нельзя таким кавалерийским наскоком решать CAS-овские задачи. Пока что у Вас получилось черт знает что, потрачена уйма машинного времени -- и все зря.

Нам нужно получить энергии S0, T1, S1 в равновесной геометрии S0 (для начала -- взятой из ХФ или DFT), те же энергии в равновесной геометрии T1 (взятой, соответственно, из ROHF или RO-DFT) и те же энергии в равновесной геометрии S1 (взятой из CIS или TDDFT). Итого -- 3 single-point CASSCF расчета с pures=.f. и усреднением по трем состояниям (S0, T1, S1). Выкладывайте выдачи и будем разбираться.

[Tonnie Fox]

А зачем Вам были нужны эти гессианы, кроме как убедиться, что у Вас минимум (или не минимум)?

Только лишь чтобы убедиться, по привычке. В Гауссиане всегда рассчитываю частоты, правда в упрощенном варианте freq=noraman, и эта процедура не занимает много времени.

Итого -- 3 single-point CASSCF расчета с pures=.f. и усреднением по трем состояниям (S0, T1, S1)

правильно ли я понял, что для этого необходимо прописать $MCSCF ntrack=3 для отслеживания трех состояний, а в $DET pures=.f. WSTATE(1)=1.0,1.0,1.0
И если требуется оптимизация, скажем, T1, то istate=2, а для S1 istate=3, при этом мультиплетность =1 в $CONTRL в обоих случаях?

[sanya1024]

Сделайте пока без оптимизации. Просто расчет энергий в точке. И тогда трекинг не нужен, и istate можно не указывать. Остальное все правильно.

[Tonnie Fox]

про оптимизацию я на будущее спросил)
сингл-поинтс рассчитываю, но приходится ставить большое значение maxit, т.к. 40-50 итерациями дело не обходится. Например, первый расчет закончился на 130+ шаге. В принципе по времени совсем немного выходит.

[sanya1024]

На будущее -- усреднять надо по тем состояниям, к-рые нас интересуют. Интересует набор S0, T1, S1 -- значит, так. Если вдруг окажется, что понадобится еще какой-нибудь триплет или синглет -- придется включить и их.

С maxit -- лучше делать не так. Лучше ставьте небольшое значение maxit и почаще рестартуйте расчет, контролируя МО (NATURAL и OPTIMIZED) на промежуточных шагах, там где MCSCF IS NOT CONVERGED. Очень может быть, что понадобится какая-то коррекция, типа выкинуть из акт. пространства МО с заселенностью почти 2 (~1.99) и/или почти 0, или наоборот, добавить орбитали.