import hfscf
import numpy as np
def SCF (r = 1.4632, Z=[1,1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3):
# Задаем центры атомов
R = [0, r]
# Генерация матрицы перекрывания S
s_scf = hfscf.S(R, b1, b2, b)
# Генерация гамильтониана H для оператора Фока (core-Hamiltonian operator)
h_scf = hfscf.H(R, Z, b1, b2, b)
# Диагонализация матрицы S и поиск диагональной матрицы X (X = S^(-1/2))
# Xa - сопряженно-транспонированная матрица
X = hfscf.diagon(m=s_scf)
Xa = X.getH()
# Создание матрицы плотности P
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
# Начинаем проверку сходимости SCF (итеративно заполняем матрицу F, диагонализируем, проверяем разницу)
for iteracion in range(50):
# Расчет двух электронной части матрицы Фока (G)
# Генерация матрицы Фока как:
# F = H + G
g_scf = hfscf.G(r, p_scf, b1, b2, b)
f_scf = h_scf + g_scf
# Создание матрицы F'
# F' = X_adj * F * X
f_tra = Xa * f_scf * X
# Диагонализация F' и генерация C'
c_tra = hfscf.diagon2(m=f_tra)
# C = X * C'
# Построение матрицы коэффициентов C
c_scf = X * c_tra
# Из матрицы C строим матрицу Р
# Пересчет матрицы плотности P
p_temp = hfscf.P(C=c_scf)
print("Завершил " + str(iteracion + 1) + " итерацию.\n")
# Проверка сходимости
if np.linalg.norm(p_temp - p_scf) < 1E-4:
print("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("SCF не сошлось.")
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
orbitals = SCF()
Завершил 1 итерацию. Завершил 2 итерацию. Завершил 3 итерацию. Завершил 4 итерацию. Завершил 5 итерацию. Завершил 6 итерацию. Завершил 7 итерацию. SCF сошлось.
orbitals
{'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 ]])}
Нарисуем орбитали в 1 измерении
hfscf.orbital(orbitals['C'])
Теперь в двух
hfscf.orbital2D(orbitals['C'], orbitals['X'], orbitals['F'])
Видны 2 орбитали (связывающая слева и разрыхляющая справа)
orben = hfscf.ener_orbs(orbitals['X'], orbitals['F'])
elecen = hfscf.ener_elec(orbitals['P'], orbitals['H'], orbitals['F'])
toten = hfscf.ener_tot()
print('Энергии орбиталей: ' + str(orben[0]) + ' и ' + str(orben[1]))
print('Электронная энергия молекулы: ' + str(elecen))
print('Полная энергия молекулы: ' + str(toten))
Энергии орбиталей: -0.5646153594583119 и 0.6368673980935621 Электронная энергия молекулы: -1.7974485548087507 Полная энергия молекулы: 0.683433570256971