Как извлечь волновые функции из результатов расчета в GAMESS

[Osmiy]

Подскажите пожалуста, где можно почитать о том, как представляются волновые функции в файлах результатов расчетов программы GAMESS (ну или дргих программ), их аналитический вид (в виде экспонент) для каждого базиса и др. Главное что бы можно было расчитать электронную плотность для молекулы.

[VTur]

Вид их таков:
ATOMIC BASIS SET
----------------
THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED
THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY

SHELL TYPE PRIM EXPONENT CONTRACTION COEFFICIENTS

C

10 S 1 4563.240000 0.778202 ( 0.001967)
10 S 2 682.024000 1.448695 ( 0.015231)
10 S 3 154.973000 2.383089 ( 0.076127)
10 S 4 44.455300 3.200089 ( 0.260801)
10 S 5 13.029000 3.013004 ( 0.616462)
10 S 6 1.827730 0.247599 ( 0.221006)

11 L 7 20.964200 0.800629 ( 0.114660) 2.573592 ( 0.040249)
11 L 8 4.803310 2.127420 ( 0.919999) 2.408248 ( 0.237594)
11 L 9 1.459330 -0.002868 ( -0.003031) 1.865281 ( 0.815854)

12 L 10 0.483456 0.413217 ( 1.000000) 0.574627 ( 1.000000)

13 L 11 0.145585 0.167976 ( 1.000000) 0.128185 ( 1.000000)

14 L 12 0.043800 0.068236 ( 1.000000) 0.028562 ( 1.000000)

15 D 13 2.504000 8.203884 ( 1.000000)

16 D 14 0.626000 0.725128 ( 1.000000)

17 D 15 0.156500 0.064093 ( 1.000000)

18 F 16 0.800000 0.891060 ( 1.000000)

Сначала стоит номер оболочки, потом тип функции (все декартовы - 6d 10f), номер гауссовой функции, показатель степени, коэффициент перед экспонентой, коэффициент перед группой (сжатых) функций.

А плотность (или, может быть, матрицу плотности?) руками считать будете?

[ИСН]

Ну, можно попросить прогу выдать плотность в файл. Поточечно.

[VTur]

Есть такая штука - cube файл, размер его меньше 20 МБ не видел. Обрабатывать руками - одно удовольствие.

[Osmiy]

А плотность (или, может быть, матрицу плотности?) руками считать будете?.

Электронную плотность буду считать сам. Сначала хочу написать программу визуализации ЭП, за одно можно правильность расчетов визуально проверить. Потом, возможно, буду писать пограмму по поиску критических точек и их параметров.

[Osmiy]

Вот сопоставляю пример из одной лекции с расчетом в базисе 3-21:


C

1 S 1 172.2560000 0.061766907377
1 S 2 25.9109000 0.358794042852
1 S 3 5.5333500 0.700713083689

2 L 4 3.6649800 -0.395895162119 0.236459946619
2 L 5 0.7705450 1.215834355681 0.860618805716

3 L 6 0.1958570 1.000000000000 1.000000000000


2px= 0.236*x*exp(-3.665r^2)+0.861*x*exp(-0.771r^2)+*"лямбда"*x*exp(-0.196r^2)


Откуда берутся цифры в общем понятно, только наверное там опечатка, и вместо *"лямбда" должно быть 1.000.
r- как я понимаю, удаленность от ядра.
х- тоже самое, что и r, только учитывается еще и знак.
Как построить py и pz я догадываюсь. А как будут выглядеть 1s и 2s?
Почему не используются числа -0.396 и 1.216?

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

И вот еще что- это только же радиальная часть функции, а угловая часть просто перется из водорода?

[Nord]

1. Прочитайте про контрактацию базисов.
2. Электронная плотность зависит не только от базиса, но и от коэффициентов, получаемых методами квантовой химии. В приведенной Вами формуле этих коэффициентов не наблюдается.
3. Угловая зависимость задается сферическими гармониками, точнее т.н. solid harmonics -- однородными многочленами от x,y и z, которые образуют некоторое неприводимое представление группы SO(3) (нумеруемое значением l оболочки). Этот результат справедлив для любой сферически-симметричной (атомной) задачи, а не только для водородоподобной.

[Osmiy]

У меня ничего нет про контрактацию базисов и в интернете ничего не нашел.
Nord, какие вы имели ввиду коэффициенты, те которые определяют степень участия АО в образовании МО.
Может вы знаете, где можно почитать нормальное руководство по данной теме.

[Nord]

Могу порекомендовать
Jensen F. Introduction to computational chemistry. Wiley, 1999, ISBN 0471980854

[Rt19]

Вот сопоставляю пример из одной лекции с расчетом в базисе 3-21:

C

1 S 1 172.2560000 0.061766907377
1 S 2 25.9109000 0.358794042852
1 S 3 5.5333500 0.700713083689

2 L 4 3.6649800 -0.395895162119 0.236459946619
2 L 5 0.7705450 1.215834355681 0.860618805716

3 L 6 0.1958570 1.000000000000 1.000000000000

2px= 0.236*x*exp(-3.665r^2)+0.861*x*exp(-0.771r^2)+*"лямбда"*x*exp(-0.196r^2)

Откуда берутся цифры в общем понятно, только наверное там опечатка, и вместо *"лямбда" должно быть 1.000.
r- как я понимаю, удаленность от ядра.
х- тоже самое, что и r, только учитывается еще и знак.
Как построить py и pz я догадываюсь. А как будут выглядеть 1s и 2s?
Почему не используются числа -0.396 и 1.216?

Выбросите лекции и убейте лектора.

Вкратце. В написанном базисном наборе задано 9 функций.
В выдаче орбиталей в out'e (хотя я рекомендую использовать $VEC в punch'e - там представление данных и разрядность выше и точнее)
будет что-то типа для каждой орбитали (a1-a9) - числа.
C S a1
C S a2
C PX a3
C PY a4
C PZ a5
C S a6
C PX a7
C PY a8
C PZ a9

Соответствующая МО имеет вид
MO=a1*(
0.061766907377*FS(172.256)*exp(-172.256*(r-rC)^2)+
0.358794042852*FS(25.9109)*exp(-25.9109*(r-rC)^2)+
0.700713083689*FS(5.53335)*exp(-5.53335*(r-rC)^2))+
a2*(
-0.395895162119*FS(3.66498)*exp(-3.66498*(r-rC)^2)+
1.215834355681*FS(0.770545)*exp(-0.770545*(r-rC)^2))
+a3*(x-xC)*(
0.236459946619*FP(3.66498)*exp(-3.66498*(r-rC)^2)+
0.860618805716*FP(0.770545)*exp(-0.770545*(r-rC)^2))+
+a4*(y-yC)*(
0.236459946619*FP(3.66498)*exp(-3.66498*(r-rC)^2)+
0.860618805716*FP(0.770545)*exp(-0.770545*(r-rC)^2))+ ... и т.д.

FS, FP и тд - соответствующие нормы примитивных функций.
int((FS(a)*exp(-a*r*r))^2) dV = 1
int((FP(a)*x*exp(-a*r*r))^2) dV = 1 и т.д.

Удачи!

P.S. Все вышенаписанное верно в случае, когда исходный базис нормирован на 1. В противном случае перед
каждым слагаемым появляется еще норма базисной функции:
MO=a1*Norm1S*(
0.061766907377*FS(172.256)*exp(-172.256*(r-rC)^2)+
0.358794042852*FS(25.9109)*exp(-25.9109*(r-rC)^2)+
0.700713083689*FS(5.53335)*exp(-5.53335*(r-rC)^2))+
a2*Norm2s*(...

int(Norm1S*(0.061766907377*FS(172.256)*exp(-172.256*(r-rC)^2)+
0.358794042852*FS(25.9109)*exp(-25.9109*(r-rC)^2)+0.700713083689*FS(5.53335)*exp(-5.53335*(r-rC)^2)))^2 dV =1

[Nord]

Вот еще хороший материал на русском языке
http://www.qchem.ru/d/khsl-qchem_lect/13Lecture-12.pdf

Upd
А, Вы, кажется, об этих лекциях и говорили...

[Osmiy]

Rt19, а как узнать нужно ли расчитывать нормы базисных функций Norm1S, Norm2S и т.д.?
И еще, как я понимаю, если для определенной МО подсчитать a1^2+a2^2+...+an^2, то должна получиться единица. Но это не получается. С этим я давно еще столкнулся.

[Nord]

И еще, как я понимаю, если для определенной МО подсчитать a1^2+a2^2+...+an^2, то должна получиться единица. Но это не получается. С этим я давно еще столкнулся.

Нет, не получится. Атомные функции на разных центрах неортогональны. Но интеграл перекрывания стремится к нулю при разведении ядер и тогда равенство выполнено. Проверьте на молекуле водорода с разведенными на 100 А ядрами.

[Rt19]

Rt19, а как узнать нужно ли расчитывать нормы базисных функций Norm1S, Norm2S и т.д.?

Вообще-то по хорошему это самому нужно сделать для всех базисных функций.
Но можно доверить это и GAMESS'у. Вот такие строчки если есть:
*** WARNING! ATOM 1 SHELL 1 TYPE S HAS NORMALIZATION 1.02152578
*** WARNING! ATOM 1 SHELL 2 TYPE S HAS NORMALIZATION 1.73798437
*** WARNING! ATOM 1 SHELL 5 TYPE P HAS NORMALIZATION 2.80994848
то они и указывают, какие конкретно функции нужно перенормировать и на какую величину.
Хотя, повторяюсь, сделать это самому ничего не стоит.

И еще, как я понимаю, если для определенной МО подсчитать a1^2+a2^2+...+an^2, то должна получиться единица. Но это не получается. С этим я давно еще столкнулся.

Nord по существу уже ответил. Правда не только для разных центрах, но и для базисных функциях на одном центре они неортогональны. Но если молекулу водорода в STO-3G развести на бесконечность, то тогда сумма действительно равна 1.

[Osmiy]

У меня последние вопросы.
1)Я взял один атом углерода и сделал для него расчет. Взял третью МО (она состоит из 2рх и 3рх) и сделал для нее визуализацию ЭП (результат в файле). Получилась нормальная р-орбиталь. Т.е. угловая часть функции уже учтена в экспонентах и сферические гармоники не нужны?
2)Как я понимаю, FS(a),FP(a) нужно расчитывать всегда?

[Rt19]

1) Вы ее уже учли, умножив на (x-xC)

x*P(r)=r*sin(theta)*cos(fi)*P1(r)=1/2(Y_{11}+Y_{1-1})*P2(r)
Вообще про угловые зависимости в виде произведения степеней x,y и z
Nord выше уже писал.

2) Зависит от программы, но де-факто стандартом - да

[Osmiy]

Я тут вот еще о чем подумал: у одного атома р-орбитали будут строго направлены вдоль осей, но в молекуле они могут быть направлены произвольно, как это нужно учитывать или это как-то уже учитывается?

[Nord]

Я тут вот еще о чем подумал: у одного атома р-орбитали будут строго направлены вдоль осей, но в молекуле они могут быть направлены произвольно, как это нужно учитывать или это как-то уже учитывается?

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

[Osmiy]

Ну что ж вроде все узнал, что нужно.
Спасибо, всем!

[VTur]

Если Вам нужны численные значения эл. плотности для визуализации и расчетов с ее использованием, то, может, действительно проще использовать cube файл. Он представляет набор записей для сетки точек: координаты точки и значение эл. плотности. В него можно включитиь проекции градиента, лапласиан, разность альфа и бета эл. плотностей, данные по орбиталям для всех точек.