Практикум 4. Вычисление параметров молекулы водорода¶
In [1]:
from hfscf import *
Расчет SCF¶
In [2]:
def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
R = [0, r]
if vbs: print("*) Построение матрицы перекрывания")
s_scf = S(R, b1, b2, b)
if vbs: print("\n*) Построение гамильтониана")
h_scf = H(R, Z, b1, b2, b)
# Диагонализация S и получение матриц X и Xa (присоединенной матрицы X)
if vbs: print("\n*) Диагонализация матрицы перекрывания")
X = diagon(m=s_scf)
Xa = X.getH()
# Создание матрицы P
if vbs: print("\n*) Построение исходной матрицы плотности")
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64) # Referencia (7) p. 148
# Начало итераций
if vbs: print("\n*) Старт итераций SCF")
for iteracion in range(50):
# Построение матрицы Фока
# F = H + G
if vbs: print("\n**) Вычисление интегралов электронов, построение матрицы Фока")
g_scf = G(r, p_scf, b1, b2, b)
f_scf = h_scf + g_scf # Referencia (7) p. 141 eq. (3.154)
# Нахождение F'
# F' = X_adj * F * X
if vbs: print("**) Построение матрицы F'")
f_tra = Xa * f_scf * X
# Диагонализация F'
if vbs: print("**) Диагонализация матрицы F' и получение C'")
c_tra = diagon2(m=f_tra)
# Расчет С
# C = X * C'
if vbs: print("**) Вычисление матрицы коэффициентов")
c_scf = X * c_tra
# Расчет P из матрицы коэффициентов C
if vbs: print("**) Пересчет матрицы плотности")
p_temp = P(C=c_scf)
print("\n Итерация " + str(iteracion + 1) + " завершена\n")
# Проверка сходимости
if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
print("\n\n--> Сходимость SCF достигнута")
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
else:
p_scf = p_temp
print("\n\n--> Нет сходимости SCF\nПересмотрите предположения")
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
In [3]:
h2_scf = SCF(vbs=True)
*) Построение матрицы перекрывания *) Построение гамильтониана *) Диагонализация матрицы перекрывания *) Построение исходной матрицы плотности *) Старт итераций SCF **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 1 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 2 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 3 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 4 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 5 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 6 завершена **) Вычисление интегралов электронов, построение матрицы Фока **) Построение матрицы F' **) Диагонализация матрицы F' и получение C' **) Вычисление матрицы коэффициентов **) Пересчет матрицы плотности Итерация 7 завершена --> Сходимость SCF достигнута
In [4]:
h2_scf
Out[4]:
{'S': matrix([[0.99999999, 0.63749012], [0.63749012, 0.99999999]]), 'H': matrix([[-1.09920375, -0.91954841], [-0.91954841, -1.09920375]]), 'X': matrix([[ 0.55258063, 1.17442445], [ 0.55258063, -1.17442445]]), 'F': matrix([[-0.34684107, -0.57771139], [-0.57771139, -0.34684028]]), 'C': matrix([[ 0.55258114, -1.17442421], [ 0.55258013, 1.17442468]]), 'P': matrix([[0.61069182, 0.61069071], [0.61069071, 0.6106896 ]])}
In [5]:
elec = ener_elec(h2_scf['P'], h2_scf['H'], h2_scf['F'])
total = ener_tot(elec=elec)
print('Электронная энергия молекулы: {:.2f} a. u. ({:.2f} ккал/моль), общая энергия: {:.2f} a. u. ({:.2f} ккал/моль)'.format(elec, elec * 627.51, total, total * 627.51))
Электронная энергия молекулы: -1.80 a. u. (-1127.92 ккал/моль), общая энергия: -1.11 a. u. (-699.06 ккал/моль)
Визуализация орбиталей¶
In [6]:
orbital(h2_scf['C'])
In [7]:
orbital2D(h2_scf['C'], h2_scf['X'], h2_scf['F'])
Итого получаем две молекулярные орбитали, которые педставляют собой линейные комбинации двух атомных. Одна из орбиталей энергетически выгодна (значение энергии отрицательно), другая — нет; это и есть т. н. связывающая и разрыхляющая МО.