расчёт комплексов органических лигандов с Pd

[VIPer]

поясните пожалуйста разницу между:
"Overlap populations:
atom charge"
и
"Atomic partitions of the density:
atom charge"
значения последних примерно в 2 раза меньше (по модулю), знак на одинаковых атомах одинаковый. понятно, что в расчётах частичные заряды не имеют того смысла, что на практике, но это какие-то разные представления частичных зарядов на атомах или вообще разные по физическому смыслу величины? какие из них тогда можно понимать, как привычные "частичные заряды на атомах" или условно считать ими?

[sanya1024]

Одни заряды (overlap populations) -- по Малликену, другие -- по Хиршфельду. Последние используются гораздо реже, поскольку не везде реализованы. Малликен -- он во всех программах есть, можно сравнивать разные расчеты. Сейчас меня тапками закидают, что, мол, у Малликеновских зарядов полным-полно недостатков, и что от базиса они зависят, и от погоды на Марсе, и от хрена с луком. Но то, что малликеновские заряды считают все программы -- перевешивает. К тому же в однотипных расчетах (один функционал, один базис) серии соединений сравнение малликеновских зарядов будет вполне корректным. Мы же не абсолютные никому не нужные цифры будем сравнивать, а тенденции в ряду.

[sanya1024]

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

Лиганд на вид вполне нормальный, жесткий -- за исключением "азотистого" хвоста (атомы 20--26, 40,41). Там может быть болтанка.

Кстати, если приведенная геометрия -- стартовая, то довольно плохая. Например, связь C-O (атомы 12-13) не копланарна связям C-C (12-4) и C-N (12-3). То же самое насчет связи N-H (атомы 25-41) -- на мой взгляд, она должна быть копланарна кольцу. Какая должна быть конфигурация связи N=N (атомы 20-26) -- я затрудняюсь даже предположить. Связи C-N (атомы 10-20 и 22-26) очевидно длинноваты.

Скорее всего это должно исправиться на первых геом. шагах.

ещё вопрос по настройкам - задал расчёт с параметрами
print=charges+bonds+vectors+molden

Вот это явно лишнее для оптимизации. Вы задали вывод орбиталей (в формате Природы и в формате Molden) на каждом геом. шагу. Конечно, это не нужно. Уберите +vectors+molden -- выдача существенно ужмется. После оптимизации можете сделать single-point расчет с +vectors+molden и посмотреть орбитали ChemCraft-ом или Molden-ом.

[sanya1024]

Да, Вам все уже и без меня хорошо рассказали

[VIPer]

да, это вроде в инпуте только после молекулярной механики оптимизация.
после природы координаты следующие:

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

   
6    -2.03850166   2.34790505   2.02573285
   6    -2.09535503   1.21669722   0.96598083
   7    -0.67235425   0.85226921   0.80566755
   6    -0.65271819   2.97208491   1.79321175
   7    -2.69560342  -0.06325281   1.43969837
   6    -1.79337165  -1.09362278   1.15856539
   6    -0.53791495  -0.55099711   0.79948070
   6    -1.96172713  -2.47495009   1.21562885
   6    -0.86068542  -3.27949233   0.91375895
   6     0.40407659  -2.73775618   0.62864100
   6     0.56717605  -1.33288752   0.54328070
   6     0.21333921   1.78389695   1.36060103
   8     1.41165764   1.63779524   1.46830829
   6    -2.75334141   1.65961405  -0.34338904
   6    -4.09073127   2.08724151  -0.33014240
   6    -4.73242901   2.44935624  -1.51451118
   6    -4.04261553   2.39496999  -2.72885455
   6    -2.71206206   1.97557238  -2.74831716
   6    -2.07025598   1.60801365  -1.56269902
   7     1.39011828  -3.70142827   0.28034484
   7     2.88113309  -2.22877582   2.44861652
   6     3.17982968  -2.45449827   1.16513689
   7     4.20729260  -1.70805300   0.65520768
   6     4.54250890  -0.94601974   1.68250443
   7     3.77527747  -1.25333643   2.75802580
   7     2.62820629  -3.52962636   0.41919759
   1    -2.08659920   1.87706529   3.02008808
   1    -2.86614407   3.06350848   1.92315597
   1    -0.21328774   3.44658544   2.68092419
   1    -0.67419000   3.71741379   0.97957927
   1    -3.65564053  -0.21134331   1.12264656
   1    -2.92610041  -2.91878516   1.47280780
   1    -0.95224158  -4.36798493   0.90422793
   1     1.51521940  -0.87366143   0.26622309
   1    -4.64006038   2.14450083   0.61573725
   1    -5.77332345   2.78166987  -1.48887227
   1    -4.54294230   2.68182824  -3.65698644
   1    -2.16519199   1.93120815  -3.69348254
   1    -1.02880498   1.27985011  -1.57785062
   1     5.31532529  -0.18006733   1.69999489
   1     3.79845115  -0.86808373   3.70081318

и вообще - как только после оптимизации запускаю с новыми координатами какой-нибудь расчёт, если слово "optimize" присутствует в параметрах, она всегда находит, что ещё "пооптимизировать", и ещё шагов 160 чего-то там выверяет))

[VIPer]

Да, Вам все уже и без меня хорошо рассказали

нет, нет, без Вас и Ваших подробных, терпеливых разъяснений не обойтись

[sanya1024]

Ага, структура уже гораздо лучше. А конфигурация N=N связи точно цис, а не транс? Может, в транс-конфигурации стабильнее будет?

[VIPer]

Ага, структура уже гораздо лучше. А конфигурация N=N связи точно цис, а не транс? Может, в транс-конфигурации стабильнее будет?

да, я тоже над этим думал.
корректно ли будет вручную задать такую конфигурацию и оптимизировать её?

[Lexx]

корректно ли будет вручную задать такую конфигурацию и оптимизировать её?

да, вполне. Потом можете сравнить, какая из них получилась устойчивее.

[amge]

а вот лиганд с такой толерантностью (-5) сошёлся через 719 шагов..

Какое-то запредельное число шагов. Что-то не так. Можно инпут? Хочется посмотреть, в чем тут дело.

VIPer, все же что-то у Вас не так. У меня геометрия из Вашего инпута (самая первая) тем же методом в том же базисе и с теми же лимитами, что были в инпуте, сошлась за 95 шагов. Конечно, и это много, но понятно, почему: фенилу пришлось на 90 гр. поворачиваться. Priroda 6 (2006.08.20).

[EvgeniX]

Странно, у меня за 250 шагов оптимизировалась.

[sanya1024]

Думаю, дело в ошибках округления (у нас, как мы помним, $grid accuracy=1e-8 стоит) в разных версиях (даже, скорее, в разных сборках) Природы на разных машинах. Я с таким сталкивалась на особо противной системе -- комплексы Eu3+ в полноэлектронном базисе. На одном кластере вообще не удавалось добиться сходимости ССП в точке, на другом сходимость (при использовании того же инпута) худо-бедно достигалась. Когда поставила $grid accuracy=1e-9, сходиться стало и там, и там -- но за разное число шагов. Потом, правда, вообще пришлось отказаться от DFT для лантанидных комплексов.

[VIPer]

у меня вот даже после прошедшей оптимизации эта структура может ещё сотни шагов оптимизироваться.. причём параметры этих структур отличаются незначительно.. такое ощущение, что алгоритм по кругу ходит практически

[sanya1024]

Grid accuracy поставили 1e-9? ставьте, глюков станет меньше. И ужесточите критерий сходимости ССП на порядок по сравнению с дефолтным. Кстати, есть еще такая полезная опция $scf iterations=n,m $end . Первая цифра n -- общее число итераций ССП, вторая m -- число т.н. микроитераций (это какие-то тонкости алгоритма, не спрашивайте у меня подробностей -- не знаю). По дефолту микроитераций делается 16, потом следующие 16 микроитераций (или достигается какой-то критерий сходимости микроитераций). Если за 16 микроитераций критерий не достигнут, а расчет продолжается, иногда это может стать источником накопления погрешностей, численных нестабильностей и болтанки -- видимо, это Ваш случай. Увеличьте число микроитераций (например, до 32) и итераций (до 100, больше смысла нет). И $scf conv=1e-7. Посмотрите, что будет.

[VIPer]

спасибо за совет!
сейчас попробую поставлю

[amge]

VIPer, сочувствую, неприятная Вам попалась модель. Видимо, действительно что-то в ней болтается настолько легко, что даже ошибки округления играют роль. Но это совершенно нетипичная ситуация. Поэтому предложение. Если у Вас нет твердой уверенности в том, что Вы имеете дело с лигандом цис-конфигурации N=N, то, может стоит перейти к лиганду с транс-конфигурацией? Последний намного (на 10 с лишним ккал/моль) стабильнее. И с ним, вдруг повезет, и не будет проблем со сходимостью. Цис-лиганд, конечно, привлекателен тем, что его азоты образуют "коробочку", в которую напрашивается встроить металл. Но и с транс-лигандом можно что-то придумать, например, металл координируется с одним местом одного лиганда и с другим местом другого лиганда, образуя, в случае соотношения металл-лиганд 1:1, что-то вроде полимера.

На 1-й картинке структура цис-лиганда, получившаяся у меня из Вашего инпута, на 2-й - несколько конформеров транс-лиганда, оптимизированные DFT/PBE/rL1, в файлах есть их энергии.

lig-c.sdf

lig-t.sdf

[sanya1024]

К тому же никто не мешает лиганду при комплексообразовании перейти в цис-форму. Для N=N связи барьер не должен быть большим.

[VIPer]

Посмотрите, что будет.

действительно, сошлось всего за 97 шагов! гораздо быстрее, безусловно.

amge, огромное Вам спасибо! обязательно буду и транс-изомер рассматривать, оптимизация его уже тоже проведена.

[VIPer]

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

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

 Electronic Spectrum:
 +-------+-------------------+------------+-------------+------------+
 | state | excitation energy | wavelength | transition  | oscillator |
 | no.   |    (au)     (ev)  |    (nm)    | moment (au) | strength   |
 +-------+-------------------+------------+-------------+------------+
 |     1 |  0.08596   2.3392 |  530.03234 |    1.039114 |   0.061880 |
 |     2 |  0.10265   2.7933 |  443.86645 |    0.091584 |   0.000574 |
 |     3 |  0.10939   2.9767 |  416.51516 |    1.249003 |   0.113768 |
 |     4 |  0.11009   2.9956 |  413.88514 |    0.286584 |   0.006028 |
 |     5 |  0.11915   3.2422 |  382.40569 |    0.397552 |   0.012554 |
 |     6 |  0.12542   3.4130 |  363.27697 |    0.079091 |   0.000523 |
 +-------+-------------------+------------+-------------+------------+

не поверите, но как раз в реальном спектре, снятом сегодня, имеется пик с максимумом где-то 412-416 нм совпадение прям точное.

посмотрим, как комплекс теперь будет поглощать

хм, и комплекс прям более, чем попал
расчёты:

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

 Electronic Spectrum:
 +-------+-------------------+------------+-------------+------------+
 | state | excitation energy | wavelength | transition  | oscillator |
 | no.   |    (au)     (ev)  |    (nm)    | moment (au) | strength   |
 +-------+-------------------+------------+-------------+------------+
 |     1 |  0.04738   1.2893 |  961.65602 |    0.114325 |   0.000413 |
 |     2 |  0.05125   1.3945 |  889.08765 |    0.351702 |   0.004226 |
 |     3 |  0.06776   1.8438 |  672.44651 |    1.116447 |   0.056305 |
 |     4 |  0.06974   1.8978 |  653.31283 |    1.515655 |   0.106808 |
 |     5 |  0.07091   1.9296 |  642.54280 |    0.675270 |   0.021556 |
 |     6 |  0.07879   2.1441 |  578.25411 |    0.228835 |   0.002751 |
 +-------+-------------------+------------+-------------+------------+

в реальном спектре его максимум лежит около 646 нм!! ну прям супер)

соотношение для величины силы осциллятора и коэффициентом экстинкции нашёл, но тогда для перевода в опт. плотность надо было конечно концентрации готовить точно, а мы их на глаз так делали..

[sanya1024]

Оптическую плотность? ишь, чего захотели Чтобы перевести силу осциллятора в коэф. экстинкции, нужна еще и ширина пика. Можно спектр, нарисованный в виде "палок" (сила осциллятора vs. длина волны) уширить Гауссом так, чтобы сила осциллятора соответствовала площади пика (как это должно быть), а не высоте (как это часто делается). При этом полуширину пика можно задавать по вкусу. Но согласие с эксп. вряд ли будет хорошее: у Вас раствор -- раз, полоса в спектре не одна (соответственно, Вы вряд ли сможете надежно отделить этот пик от соседнего и точно определить его площадь) -- два, ну, и наверное я что-нибудь еще забыла.
Идеальное совпадение посчитанных в газовой фазе спектров с экспериментальными значениями для раствора -- повод насторожиться У меня, например, первая мысль была: Вам как новичку повезло. Вторая: такое везение тянет на лишние вопросы на защите, особенно если попадется скептик, не доверяющий расчетам в принципе. Лучше не представлять этот результат как великое достижение. Кстати, какой был функционал и в каком растворителе снимали спектры?
Еще неплохо посмотреть орбитальный состав интересующего перехода. Это запись типа

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

 state   1: E =  0.07436597

 160 a ->  161 a :  0.9884031
 160 a ->  162 a :  0.1318887 

Чтобы посмотреть на собс-но орбитали, нужно сделать single-point расчет (task=energy), но написать print=+vectors+molden, потом вырезать из полученной здоровенной выдачи все строчки, начинающиеся на mos>, а потом это слово mos> из всех строчек удалить -- это и будет молденовская выдача -- дайте ей расширение *.molden. Ее можно посмотреть ChemCraft-ом.