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

[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 $$

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