Решение ядерного уравнения Шредингера (ЖРГО)

[B3LYP]

Мне нужно для своих задач научиться решать ядерное уравнение Шредингера в приближении “Жёсткий ротатор – гармонический осциллятор”. Изначально для этого надо найти поверхность потенциальной энергии молекулы в гармоническом приближении. Я провёл такой расчёт для молекулы озона в программе Gaussian09. Программа мне выдала выходной файл (с расширением .out) и файл “форматированный чекпоинт” с расширением .fch.

В файле .fch напечатаны такие декартовы координаты атомов:

O1 0.000000000 0.000000000 0.455483035
O2 0.000000000 1.131064087 -0.227741017
O3 0.000000000 -1.131064087 -0.227741017

Далее приводятся силы (градиенты энергии):

Цитата:

Cartesian Gradient R N= 9
1.35706545E-15 -6.34041430E-16 9.06506923E-08 -1.71369282E-16 7.92786370E-08
-4.53253417E-08 -1.18569617E-15 -7.92786725E-08 -4.53253275E-08

Я привёл этот кусок чтобы проиллюстрировать, что моя молекула находится в положении минимума энергии (равновесная геометрия) – градиенты примерно нулевые.
Далее в файле напечатаны силовые постоянные (вторые производные энергии по декартовым координатам атомов):

Цитата:

Cartesian Force Constants R N= 45
-1.26603926E-04 3.64626364E-12 3.53287483E-01 -3.87525061E-12 -1.25555859E-12
3.67103438E-01 6.33019612E-05 3.70967163E-12 3.02013391E-12 -7.00734661E-05
-5.91922853E-12 -1.76643742E-01 1.18068182E-01 1.88796760E-13 2.81060001E-01
3.22031734E-12 1.06807784E-01 -1.83551719E-01 -3.06840880E-12 -1.12437983E-01
1.23835123E-01 6.33019604E-05 6.01812544E-12 -1.98130066E-12 6.77151024E-06
-3.91208118E-12 2.88949326E-13 -7.00734673E-05 -5.40095541E-15 -1.76643742E-01
-1.18068182E-01 -3.91267352E-12 -1.04416259E-01 5.63019923E-03 1.86526736E-13
2.81060001E-01 -3.82414369E-13 -1.06807784E-01 -1.83551719E-01 -2.89212401E-13
-5.63019923E-03 5.97165964E-02 3.06718627E-12 1.12437983E-01 1.23835123E-01

Чтобы было понятнее, матрица вторых производных для этой молекулы выглядит так:

-0.0001 0.0000 0.0000 0.0001 0.0000 0.0000 0.0001 0.0000 0.0000
0.0000 0.3533 0.0000 0.0000 -0.1766 0.1068 0.0000 -0.1766 -0.1068
0.0000 0.0000 0.3671 0.0000 0.1181 -0.1836 0.0000 -0.1181 -0.1836
0.0001 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 -0.1766 0.1181 0.0000 0.2811 -0.1124 0.0000 -0.1044 -0.0056
0.0000 0.1068 -0.1836 0.0000 -0.1124 0.1238 0.0000 0.0056 0.0597
0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000
0.0000 -0.1766 -0.1181 0.0000 -0.1044 0.0056 0.0000 0.2811 0.1124
0.0000 -0.1068 -0.1836 0.0000 -0.0056 0.0597 0.0000 0.1124 0.1238

Такая же матрица показана в .out файле, так что вроде всё правильно.
Как мне объяснили, матрицу силовых постоянных надо диагонализировать, тогда диагональными элементами станут фундаментальные частоты колебаний данной молекулы (точнее, квадраты циклических частот). Я продиагонализировал и к сожалению не вышло ничего путного. Вот какие элементы на диагонали получилось после диагонализации:

-0.0001 0.3533 0.3671 -0.0001 0.2811 0.1238 -0.0001 0.2811 0.1238

Тут явно вышла ерунда поскольку ненулевые значения имеют разный знак. Согласно расчёту в программе Gaussian, эта молекула должна иметь такие частоты колебаний: 639.4649, 991.0042, 1066.5987.
Вначале мне сказали что моя ошибка в том, что я не переводил матрицу силовых постоянных в масс-взвешенные координаты; но когда начал переводить – лучше не стало, а потом я взял именно молекулу озона, у которой все атомы одинаковы, соответственно с точностью до множителя всё должно быть корректно и без этого перевода координат.
Ещё я возможно слышал, что для этой задачи надо ещё использовать G-матрицу, но я пока не понимаю что это такое.
Прошу помочь.

[madschumacher]

B3LYP в сообщении #1644851 писал(а):

Прошу помочь.

Я же вроде уже писал о том, как надо это считать на кемпорте? Ну раз уж здесь больше возможностей для размещения формул и кода, то вот как оно работает.

• Сначала нужно взять матрицу силовых постоянных $\mathcal{F}$ размера $3N \times 3N$, где N -- число атомов (в данном случае это N=3). Это Вы сделали правильно, и матрица приведённая выше -- правильная. Эта матрица силовых постоянных дана в атомных единицах, т.е. элементы этой матрицы $F_{ij}$ ($i,j=1,2,\ldots, 3N$) имеют размерность $\text{хартри}/\text{бор}^2$.
• Дальше эту матрицу переводим в масс-взвешенные координаты, получая новую матрицу $\mathcal{F}^{\mathrm{MW}}$ с элементами
  
  $$ F_{ij}^{\mathrm{MW}}= \frac{F_{ij}}{\sqrt{m_i \cdot m_j}} $$
  
  Здесь $m_i$ и $m_j$ соответствуют массе атомов соответствующих координатам номер i и j. Номер конкретного атома $n=1,2,\ldots, N$ получить очень просто, достаточно поделить индексы i и j на три и взять целую часть (например, $n = \lfloor i/3 \rfloor$).
• После этого надо диагонализовать матрицу $F_{ij}^{\mathrm{MW}}$ решив задачу на собственные вектора $\mathbf{u}$ и собственные значения $\lambda$, т.е.
  
  $$ F_{ij}^{\mathrm{MW}} \mathbf{u}_k = \lambda_k \mathbf{u}_k $$
  
  Пара собственного вектора $\mathbf{u}_k $ и собственного значения $\lambda_k$ ($k=1,2,\ldots, 3N$) даёт форму нормального колебания и значение угловой частоты соответствующего колебания, т.е.
  
  $$\lambda_k = \omega_k^2$$
• Ну а дальше остаётся только перевести значение в правильные единицы (обратные сантиметры), в частоты $\tilde{\nu}$. Для этого нужно их посчитать как
  
  $$ \tilde{\nu}_k = \frac{\sqrt{\lambda_k}}{2 \cdot \pi \cdot c \cdot \tau_\mathrm{au}} = \frac{{\omega_k}}{2 \cdot \pi \cdot c \cdot \tau_\mathrm{au}} \approx 219475 \times \omega_k $$
  
  где $c=29979245800$ это скорость света в сантиметрах в секунду, а $\tau_\mathrm{au}=2.4188843265864 \times 10^{-17}$ секунд -- это единица времени в атомных единицах.

[madschumacher]

B3LYP в сообщении #1644851 писал(а):

Чтобы было понятнее, матрица вторых производных для этой молекулы выглядит так:

Чтобы не быть голословным, приведу также пример кода на Питоне, как получить из этого результата конфетку. Возьмём приведённую матрицу и сохраним её в виде файла h2o.ff.

(содержание файла h2o.ff)

-0.0001 0.0000 0.0000 0.0001 0.0000 0.0000 0.0001 0.0000 0.0000
0.0000 0.3533 0.0000 0.0000 -0.1766 0.1068 0.0000 -0.1766 -0.1068
0.0000 0.0000 0.3671 0.0000 0.1181 -0.1836 0.0000 -0.1181 -0.1836
0.0001 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 -0.1766 0.1181 0.0000 0.2811 -0.1124 0.0000 -0.1044 -0.0056
0.0000 0.1068 -0.1836 0.0000 -0.1124 0.1238 0.0000 0.0056 0.0597
0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000
0.0000 -0.1766 -0.1181 0.0000 -0.1044 0.0056 0.0000 0.2811 0.1124
0.0000 -0.1068 -0.1836 0.0000 -0.0056 0.0597 0.0000 0.1124 0.1238

Дальше просто пишем следующий скрипт:

Используется синтаксис Python

import numpy as np

# step 1: load force-field matrix
ff = np.loadtxt("h2o.ff")

# step 2.1: define vector of masses with length of 3*N (N is the number of atoms, N=3 in this case)
# atom 1 is oxygen, atoms 2 and 3 are the hydrogens
# masses of isotopes are taken from https://physics.nist.gov/cgi-bin/Compos ... d_alone.pl
mO = 15.99491461957
mH = 1.00782503223
#m = np.array([mO, mO, mO, mH, mH, mH, mH, mH, mH]) # this is for water
m = np.array([mO, mO, mO, mO, mO, mO, mO, mO, mO])  # this is for ozone

# step 2.2: convert atomic masses from atomic mass units (daltons) into atomic units
#           by multiplying with the mass of the electron
m *= 1822.888486209

# convert the force-field into the mass-weighted coordinates by dividing each ij-th element with sqrt(mi*mj),
# where mi and mj are the masses corresponding to these coordinates
for i in range(0,len(m)):
    for j in range(0,len(m)):
        ff[i][j] /= np.sqrt(m[i]*m[j])
        
# step 3: find the eigenvalues of the force-field matrix in mass-weighted coordinates
# this will give the squared angular frequencies of the vibrations in atomic units
w2 = np.linalg.eigvalsh(ff)

# step 4.1: convert these values into the normal frequencies in Hz:
timeinau = 2.4188843265864e-17 # this is the unit of time in atomic units in seconds
c = 29979245800 # this is the speed of light in cm/s
Const = 1./(2.*np.pi*timeinau*c) # and this is the conversion constant
#print(np.round(Const))
v = np.sqrt(np.array(w2, dtype=complex))* Const

# now, print the frequencies
for i,vi in enumerate(v):
    print(" %2i : %7.2f + i %7.2f cm^{-1}" % (i+1,np.real(vi),np.imag(vi)))
 

При его запуске в командной строке (например как python <имя скрипта>) мы получаем следующую выдачу:

  1 :    0.00 + i   26.47 cm^{-1}
  2 :    0.00 + i   19.97 cm^{-1}
  3 :    0.00 + i   12.85 cm^{-1}
  4 :    0.00 + i   12.85 cm^{-1}
  5 :    8.27 + i    0.00 cm^{-1}
  6 :   12.85 + i    0.00 cm^{-1}
  7 :  639.51 + i    0.00 cm^{-1}
  8 :  990.97 + i    0.00 cm^{-1}
  9 : 1066.63 + i    0.00 cm^{-1}

Первые шесть самых маленьких (и частично мнимых) частот -- это трансляция и вращение молекулы озона, а последние три -- это как раз нужные цифры.

Вот ещё.

B3LYP в сообщении #1644851 писал(а):

Вот какие элементы на диагонали получилось после диагонализации:
-0.0001 0.3533 0.3671 -0.0001 0.2811 0.1238 -0.0001 0.2811 0.1238

Тут с диагонализацией явно проблемы. Собственные значения той матрицы, что Вы привели выше следующие (посчитаны тем же NumPy-ем):

[-4.24158400e-04 -2.41421356e-04 -1.00021528e-04 -1.00000000e-04
  4.14213562e-05  1.00000000e-04  2.47550765e-01  5.94424158e-01
  6.88649257e-01]

Отсюда вполне себе уже можно вытащить частоты. Достаточно поделить это всё хозяйство на $16 \times 1823$ (масса кислорода в атомных единицах), взять корень из получившегося числа, после чего домножить его на коэффициент перевода в обратные сантиметры (219475, см. выше). На примере самой больших из значений. Для самой большой частоты (6.88649257e-01=0.688649257) мы получим:

$$ \tilde{\nu} = 218475 \times \sqrt{\frac{0.688649257}{16 \times 1823}} = 1061.6 \ [\text{см}^2] $$

Для средней (5.94424158e-01=0.594424158) результатом будет

$$ \tilde{\nu} = 218475 \times \sqrt{\frac{0.594424158}{16 \times 1823}} = 986.2 \ [\text{см}^2] $$

Для низшего по частоте же колебания (2.47550765e-01=0.247550765) мы получим

$$ \tilde{\nu} = 218475 \times \sqrt{\frac{0.247550765}{16 \times 1823}} = 636.5 \ [\text{см}^2] $$

[madschumacher]

madschumacher в сообщении #1645699 писал(а):

B3LYP в сообщении #1644851 писал(а):

Прошу помочь.

После этого надо диагонализовать матрицу $F_{ij}^{\mathrm{MW}}$ решив задачу на собственные вектора $\mathbf{u}$ и собственные значения $\lambda$, т.е.

$$ F_{ij}^{\mathrm{MW}} \mathbf{u}_k = \lambda_k \mathbf{u}_k $$

Нашёл ошибку в оформлении ( ), которую не могу исправить ( ). Матрица $\mathcal{F}^{\mathrm{MW}}$ состоит из элементов $F_{ij}^{\mathrm{MW}}$, поэтмоу задача на собственные вектора/значения имеет вид

$$ \mathcal{F}^{\mathrm{MW}}\mathbf{u}_k = \lambda_k \mathbf{u}_k $$

В остальном всё правильно.

[B3LYP]

topic158289.html

[B3LYP]

Я осилил диагонализацию матрицы (реализовал диагонализацию Якоби), и сейчас у меня такой вопрос: в каких единицах после диагонализации получаются собственные вектора данной матрицы (гессиана). По-моему они близки к нормальным модам колебаний, но всё-таки не совпадают с ними, их наверно надо ещё как-то пересчитать, может разделить на корни из масс атомов, чтобы получились нормальные моды?

[Cos(x-pi/2)]

B3LYP в сообщении #1648968 писал(а):

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

Компоненты собственных векторов это безразмерные числа. Для каждого собственного вектора они подчиняются условию нормировки собственного вектора на единицу (безразмерную). Стандартные программы, вычисляющие собственные векторы, выдают их уже таким образом нормированными.

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

[B3LYP]

Извиняюсь что по любому вопросу прошу консультацию знатоков, так уж я привык. Сейчас у меня такие вопросы. Мне раньше казалось, что программа Gaussian, после того как диагонализирует гессиан, далее ещё избавляется от поступательных и вращательных частот. Вот фрагмент выдачи Gaussian-а для молекулы воды:

Цитата:

Full mass-weighted force constant matrix:
Low frequencies --- -31.5792 -21.3113 -17.1062 -0.0017 -0.0002 -0.0002
Low frequencies --- 1619.0486 3616.3668 3781.4919

Тут поступательные или вращательные частоты по величине доходят аж до 30 обратных сантиметров, что вообще весьма много. Вроде, если не путаю, в мануале написано что эти частоты (поступательные и вращательные) определяются сравнением с тензором момента инерции молекулы.
Но далее мне непонятно: в выдаче этого расчёта далее приводятся “правильные” частоты колебаний, без вращательных и поступательных:

Цитата:

Harmonic frequencies (cm**-1), IR intensities (KM/Mole), Raman scattering
activities (A**4/AMU), depolarization ratios for plane and unpolarized
incident light, reduced masses (AMU), force constants (mDyne/A),
and normal coordinates:
1 2 3
A1 A1 B2
Frequencies -- 1619.0486 3616.3668 3781.4919
Red. masses -- 1.0918 1.0368 1.0856
Frc consts -- 1.6862 7.9893 9.1464
IR Inten -- 68.4059 4.4038 8.8273
Atom AN X Y Z X Y Z X Y Z
1 8 0.00 0.00 0.07 0.00 0.00 0.04 0.00 0.07 0.00
2 1 0.00 -0.38 -0.59 0.00 0.61 -0.35 0.00 -0.57 0.41
3 1 0.00 0.38 -0.59 0.00 -0.61 -0.35 0.00 -0.57 -0.41

Здесь мне непонятно то, что эти частоты колебаний полностью совпадают с предыдущей выдачей; а я раньше думал что после обнуления поступательных и вращательных компонент они чуть поменяются. Или я всё путаю?

[B3LYP]

Я продвигаюсь с этой задачей, но снова возникают вопросы, и хочется помощи.
Как обычно я обращаюсь и к физикам, и к математикам, и к химикам-квантовикам. От физиков мне хочется, чтобы помогли понять, как решается УШ для многомерного гармонического осциллятора, т.е. молекулы. Приближение ЖРГО подразумевает, что решается одновременно задача колебательного движения ядер в молекуле, и вращательно-поступательного.
В предыдущем сообщении я привёл не совсем удачный пример, на самом деле частоты должны немного отличаться в выдаче: сначала программа просто диагонализует гессиан и получает частоты, а дальше она находит шесть поступательных и вращательных частот, и делает их нулевыми, при этом остальные частоты (колебательные) чуть уточняются.
Мне надо понять, как это делается. Информация приведена здесь:

https://gaussian.com/vib/

После первой диагонализации, надо определить центр координат молекулы и тензор момента инерции, т.е. найти систему координат по осям вращения. Далее считается некая матрица D; первые 6 столбцов этой матрицы это нормальные моды поступательных и вращательных степеней свободы. Далее проводится какая-то ортогонализация Шмидта; если я правильно понял, находятся остальные 3N-6 столбцов, они могут быть любыми, лишь быть ортогональными первым шести столбцам. Я правильно понял?
И вроде далее уже просто – Гессиан умножается на матрицу D, точнее так:

$F_{int}=D^t F D$

Тут вроде как Гессиан переводится в так называемые внутренние координаты, хотя может я что-то пока не понял. Правильно ли я понимаю, что Гессиан во внутренних координатах будет диагонализован, и на первых шести ячейках диагонали будут нули?
Вопрос по общей математической теории этих преобразований: правильно ли я понял, что Гессиан можно представить как 3N векторов размерностью 3N каждый, и диагонализация поворачивает эти 3N векторов, т.е. умножает их на “матрицу вращения” размерностью $3N \cdot 3N$?
Дальше у меня будет такая задача - как восстановить Гессиан, если известны частоты и нормальные моды колебаний, но без 6 поступательных и вращательных компонентов.

[Утундрий]

B3LYP
Вот вы сейчас зарабатываете. А как только поймёте, так сразу перестанете зарабатывать. Оно вам надо?