import hfscf
from hfscf import *
Функция расчёта орбиталей на основе подхода Рутана Холла включает: S - матрицу перекрывания; F - матрицу Фока; X - диагонализированную матрица S; H - гамильтониан ядра; C - матрицу коэффициентов.
def SCF (r = 1.4632, Z=[1,1], b1 = GTO["H"], b2 = GTO["H"], b = 3, vbs=False):
#степени свободы, r – электронные, R – ядерные
R = [0, r]
if vbs: print("Генерируем матрицу перекрываний S.")
s_scf = S(R, b1, b2, b)
if vbs: print("Генерация гамильтониана H.")
h_scf = H(R, Z, b1, b2, b)
if vbs: print("Диагонализация матрицы S и поиск диагональной матрицы X.")
X = diagon(m=s_scf)
Xa = X.getH()
if vbs: print("Создаем матрицу плотности P.")
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64) # Referencia (7) p. 148
#запуск итеративных процессов
if vbs: print("Начинаем считать SCF.")
for iteracion in range(50):
#строим матрицу Фока F
# F = H + G
if vbs: print("Генерация матрицы Фока: вычисление интегралов двух электронов.")
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
if vbs: print("Диагонализация F' и генерация C'.")
c_tra = diagon2(m=f_tra)
#C = X * C'
if vbs: print("Построение матрицы коэффициентов C.")
c_scf = X * c_tra
#строим матрицу P на основе матрицы C
if vbs: print("Расчет матрицы плотности P.")
p_temp = P(C=c_scf)
#построим графики
orbital(c_tra, r, b1, b2, b)
orbital2D(c_scf, X, f_scf, r, b1, b2, b)
print("\n Итерация " + str(iteracion + 1) + ". завершена.\n")
#проверяем сходимость
if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
print("\n\n-->Нашлось самосогласованное поле!")
elecE = ener_elec(p_temp, h_scf, f_scf)
totE = ener_tot(r, Z, elecE)
print("Электронная энерегия:", elecE)
print("Полная энергия:", totE)
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-->Самосогласованное поле НЕ нашлось!\n Изменить допущения.")
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
SCF()