Оптимизация комплекса Ni

[AlexKsen]

Уважаемые форумчане, пытаюсь оптимизировать комплекс Ni. Сначала для первичной оптимизации пытался DFT/B3LYP/3-21G, не получилось, потом хотя бы HF/3-21G, опять ничего! (Хотя PM6 оптимизировалось). Помогите разобраться в чем причина. Спасибо!

[endlesslake]

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

[sanya1024]

Совершенно верно. В PM6 и др. подобной полуэмпирике используется нечто вроде усреднения по состояниям (не столь корректное, как в MCSCF, конечно), поэтому там все получается.
AlexKsen, у Вас с самого начала много ошибок в инпуте. (1) дефолтное кол-во итераций ССП (maxit в $contrl) -- 30, этого явно недостаточно для такого монстра даже при хорошей сходимости. (2) НАСТОЯТЕЛЬНО РЕКОМЕНДУЕТСЯ оптимизировать во внутренних координатах (посмотрите форум, я неоднократно приводила соответствующие строчки), а не в декартовых по дефолту. С остальным будем разбираться, когда возьмете нормальный базис, в минимальном базисе ловить нечего.
Сходимость ССП у Вас вообще ни к черту с первых же итераций, сейчас увеличивать maxit -- только зря гонять машину. Действительно, дело в том, что у Вас неправильная мультиплетность. Скорее всего, комплекс высокоспиновый, и понадобится ROHF (не UHF!).
И еще совет: сначала разрежьте комплекс пополам (удалите метиленовый мостик, поставьте вместо него H) и разберитесь с одной половинкой.

[AlexKsen]

Попробовал оптимизировать с мультиплетностью 5, правда базис оставил 3 - 21G, взял ROHF, но число итераций не изменял и вид координат (за это прошу прощения!). Подумал, что может получится. В итоге нет. Прикрепляю inp и out. Посмотрите, пожалуйста. А также прикрепляю файл inpnew, где я беру уникальные координаты, меняю число итераций, если Вас не затруднит, посмотрите пожалуйста и скажите все ли так.
Спасибо!

[sanya1024]

Ну, что сказать?
30 ССП итераций мало для такого монстра, ставьте 100. Если и 100 будет недостаточно, тогда действительно лучше не упираться, а посмотреть на управление сходимостью (группа $SCF).
Стартовые орбитали, видимо, вообще ни к черту. Возможно, придется попробовать стартовать с орбиталей ближайшей системы с замкнутой или наоборот, высокоспиновой открытой оболочкой.
От базиса 3-21 надо отказаться. Во-первых, Вы все равно не будете с ним работать, перейдете на 6-31 как минимум. Во-вторых, даже если ценой неимоверных усилий Вы сведете ССП в 3-21, Вы не сможете стартовать расчет в бОльшем базисе с этих МО -- не предусмотрено такого в Гамессе. В третьих, из-за своей негибкости этот базис может просто-напросто дать неправльные МО (возможно, отсюда и несходимость).
Теперь по инпуту. Новый от старого отличается только добавкой COORD=UNIQUE -- но это и так по дефолту, и DFTTYP=B3LYP -- а в этом пока смысла нет все равно.

Итак: увеличить maxit до 100 в $CONTRL, сменить базис на 6-31 (+ LANL ECP для металла). Или, если совсем лень разбираться в задании ECP, просто поставьте PP=SBKJC в $CONTRL и GBASIS=SBKJC в $BASIS. По крайней мере, валентные оболочки будут прилично описываться.
Если будут продолжаться проблемы со сходимостью (скорее всего будут), будем искать хорошее стартовое приближение. RUNTYP пока поставьте ENERGY, до оптимизации геометрии нам еще как до Шанхая. Когда ССП сведем, тогда и скажу, как задавать делокализованные координаты для оптимизации (или гляньте по форуму, неоднократно говорила).

Для управления сходимостью добавьте в $SCF DIIS=.t. DAMP=.t. SHIFT=.t. Посмотрим, что с этим получится.

ПОЖАЛУЙСТА, СДЕЛАЙТЕ ВСЕ СРАЗУ: БАЗИС, MAXIT И КЛЮЧИ В ГРУППЕ $SCF. По одной опции -- не сработает.

[AlexKsen]

Спасибо! Попробую сделать, как Вы советуете.

[AlexKsen]

Здравствуйте, сделал maxit 100, базис 6-31G для начала и поставил ключи в SCF. Кстати, по ключам у меня тоже вопрос:) DIIS=.T. - это опция полезна для сложно сходящихся случаев и отвечает за переключение между SOSCF <-> FULLNR, а для чего нужны DAMP=.t. SHIFT=.t. ? Псевдопотенциалы пока не использовал. Если использовать нужно таким образом прописывать группу ECP: $ECP Ni-ECP SBKJC H-ECP NONE O-ECP NONE C-ECP NONE $END или для всех атомов в молекуле, а как быть с базисом, т.е. вместо к примеру 6-31G, нужно взять GBASIS=SBKJC Спасибо!

[sanya1024]

Здравствуйте, сделал maxit 100,

Где? нет этого в Ваших файлах. Как было 30 итераций, так и осталось. Надо поставить maxit=100 в $CONTRL.

Правда, сейчас это бесполезно. Все равно энергия у Вас вверх ползет. Может, попробуйте $guess guess=hcore $end? А еще лучше -- попробуйте стартовать с МО ближайшей системы с замкнутой оболочкой (например, ободрав с обоих атомов Ni все неспаренные электроны и поставив соответствующий заряд (+4). Да, когда будете это делать, поставтье RUNTYP=ENERGY, молекулу с таким зарядом оптимизировать НЕ НУЖНО! Когда получите стартовые орбитали для такой (искусственной) системы, скопируйте группу $VEC в инпут для молекулы с нормальным зарядом и поставьте $guess guess=moread $end. И не забудьте про maxit=100 в $CONTRL.

базис 6-31G для начала и поставил ключи в SCF. Кстати, по ключам у меня тоже вопрос:) DIIS=.T. - это опция полезна для сложно сходящихся случаев и отвечает за переключение между SOSCF <-> FULLNR,

И тут нет. DIIS и SOSCF -- это альтернативные методы достижения сходимости. Т.е., или DIIS, или SOSCF. Где-то Вы не туда посмотрели в мануале... FULLNR вообще не используется в обычном SCF, а только в MCSCF (до к-рого мы с Вами еще не скоро дойдем), да и то ею там никогда не пользуются.

а для чего нужны DAMP=.t. SHIFT=.t. ?

Не поленитесь заглянуть в мануал, описание группы $SCF. Ну что я Вам буду все это пересказывать?

Псевдопотенциалы пока не использовал. Если использовать нужно таким образом прописывать группу ECP: $ECP Ni-ECP SBKJC H-ECP NONE O-ECP NONE C-ECP NONE $END или для всех атомов в молекуле, а как быть с базисом, т.е. вместо к примеру 6-31G, нужно взять GBASIS=SBKJC Спасибо!

Я советовала SBKJC исключительно из тех соображений, что его можно задать одним ключом, точнее, двумя. Если все-таки хотите 6-31 на всем, кроме никеля, то для никеля возьмите LANL ECP и соответствующий DZ базис (полностью это называется LANL2DZ). В этом случае группа $BASIS не нужна, но базис надо указывать для каждого атома персонально, в группе $DATA. Примерно так:

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

EUROPIUM   63.0      0.048561167        -0.045387879        -0.838956447
 S   3
  1    0.7363006000D+01 -0.1893069969D+00
  2    0.5973380000D+01  0.5419949911D+00
  3    0.2492572000D+01 -0.1305799979D+01
 S   1
  1    0.6217810000D+00  0.1000000000D+01
 S   1
  1    0.3050070000D+00  0.1000000000D+01
 S   1
  1    0.5744400000D-01  0.1000000000D+01
 S   1
  1    0.2669600000D-01  0.1000000000D+01
 P   3
  1    0.4844781000D+01  0.1552890048D+00
  2    0.2980927000D+01 -0.5288870165D+00
  3    0.7304400000D+00  0.1192228037D+01
 P   1
  1    0.3107740000D+00  0.1000000000D+01
 P   1
  1    0.9628600000D-01  0.1000000000D+01
 P   1
  1    0.3375300000D-01  0.1000000000D+01
 D   3
  1    0.2534693000D+01 -0.5829699802D-01
  2    0.7843880000D+00  0.3814979870D+00
  3    0.3320980000D+00  0.7067249760D+00
 D   1
  1    0.1307590000D+00  0.1000000000D+01
 D   1
  1    0.4939000000D-01  0.1000000000D+01

SULFUR     16.0      1.531111946         3.709648519         2.282757942
 N31        6
 D       1
    1         0.6500000000  1.00000000

Обязательно после базиса поставьте пустую строку.

Итак: (1) псевдопотенциал (SBKJC для всех или LANL2DZ для Ni и 6-31 для всех остальных) + MAXIT=100 + ICHARG=4 MULT=1 SCFTYP=RHF.
Затем (2) с полученными орбиталями ICHARG=0 MULT=5 SCFTYP=ROHF GUESS=MOREAD + $VEC в конце инпута + MAXIT=100

И смотрите не перепутайте

[AlexKsen]

sanya1024, я сделал, как Вы советовали. Получил стартовые орбитали. У меня очередной вопрос (Замучил, я Вас уже:() Нужно ли в $guess guess=moread $end прописывать NORB? Это значение будет равно 32 для моего случая. Спасибо Вам большое, что никогда не отказываете в помощи и ценном совете!

[sanya1024]

Дефолтное кол-во считываемых орбиталей (NORB) равно числу занятых орбиталей в ССП (у Вас это далеко не 32) или числу остовных+активных в CASSCF/QDPT. Во многих случаях лучше считывать все орбитали, а не только занятые. Тогда NORB надо ставить равным полному числу орбиталей.

[sanya1024]

А где там псевдопотенциал? все-таки решили мучиться с 6-31G на всех атомах? за что ж Вы себя так? да и машинку жалко... Впрочем, теперь раз уж все начали в 6-31 делать, так и продолжайте, орбитали, полученные в одном базисе, перетаскивать в другой нельзя. Правда, 6-31 без поляризационных функций (ndfunc и npfunc) довольно плох даже с качественной точки зрения. Может быть, проще повторить генерацию стартовых орбиталей в базисе 6-31G(d,p)?

И где оптимизация в DLC? у Вас оптимизация будет долго и безутешно болтаться в декартовых координатах. Результата не дождетесь никогда.

Скажите честно, какой программой Вы генерите инпут? может, лучше это руками делать, в текстовом редакторе? только группу $DATA вставляйте из КемКрафта.

[AlexKsen]

Я еще раз сделаю генерацию стартовых орбиталей, только в 6-31G(d,p). Чтобы включить оптимизацию в DLC нужны ключевые слова $zmat dlc=.t. auto=.t.? Инпут я делаю сам в текстовом редакторе.

[sanya1024]

Ага, все правильно. Только для генерации DLC еще обязательно должен в $CONTRL присутствовать ключ nzvar=<любое ненулевое число>

[АлександрI]

Гм, что это за комплекс такой, что с ним столько проблем...
Хотелось бы взглянуть не него в классическом представлении - ну там символы атомов, палочки-чёрточки или ещё как-то? На кемкрафте у меня не получается разглядеть его: только атомы и пара связей...
Может быть, исходная структура комплекса задана неправильно, откуда вы её взяли, или она умозрительна из общих соображений?

[endlesslake]

да не, там вроде надо только в кемкрафте rebond нажать

[АлександрI]

Спасибо, endlesslake.
Далее, я удалил все метильные группы и вот, он, уже не монстр (как отметила Sanya), а красавчик, с порфириновыми циклами...
Вместо метилов надо поставить водороды и оптимизация может быть более успешной.
А затем понемногу добавлять эти самые метилы...
Тогда можно будет обойтись без разрыва димера удалением метиленовых групп и оптимизировать любой из кусков, а потом их снова связать - это был бы самый чистый вариант (с метилами или без).
Правда я не совсем понял, зачем при 528 электронах в исходном комплексе (взято из out.txt топикстартера) обязательно брать мультиплетность 5 и ROHF, а не по дефолту и RHF... Ведь химические и спектральные свойства десятков подобных соединений уже описаны, начиная с хлорофила...

[sanya1024]

Ну, первое подозрение при виде переходного металла -- "а вдруг он высокоспиновый". Это надо проверить, тем более что RHF и mult=1 после мучений с ROHF и mult=5 покажется уже совсем-совсем нестрашным. Это, конечно, не порфириновый цикл, координация у Ni -- тетраэдр (что там нам говорит теория поля лигандов по этому поводу?)

[AlexKsen]

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

[sanya1024]

Ну да. Только очень прошу, архивируйте выдачи, 31 Мб -- это тяжело тащить.
Гм, похоже, с высокоспиновостью мы переборщили , кривовато выглядит сам комплекс, и граничные орбитали на один бок перекошены. Самое главное -- на металле ничего нет, это подозрительно (базис неадекватный?). Наверное, стоит попробовать на Ni дать LANL2DZ c ECP и все переделать.