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_scf = H(R=R, Z=Z, b1=b1, b2=b2, b=b)
if vbs: print("Диагонализируем матрицу S и получаем транспонированную матрицу")
X = diagon(m=s_scf)
Xa = X.getH()
if vbs: print("Считаем матрицу расчетной плотности P")
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
# улучшаем матрицу P в каждом новом цикле
if vbs: print("Подсчет SCF")
for iteration in range(50):
# F = H + G
if vbs: print("Генерируем матрицу Фока, вычисляем интеграл для двух электронов")
g_scf = G(r=r, p=p_scf, b1=b1, b2=b2, b=b)
f_scf = h_scf + g_scf
# генерируем матрицу F'
# F' = X_adj * F * X
if vbs: print("Смена базиса матрицы F")
f_tra = Xa * f_scf * X
# диагонализируем построенную матрицу F' и генерируем матрицу C'
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.")
p_temp = hfscf.P(c_scf)
#строим графики
orbital(c_tra, r, b1, b2, b)
orbital2D(c_scf, X, f_scf, r, b1, b2, b)
print("Цикл " + str(iteration + 1) + " закончен.")
# проверяем сходимость
if np.linalg.norm(p_temp - p_scf) < 1E-4:
print("\n\n-->Самосогласованное поле сходится!")
print({"S" : s_scf, "H" : h_scf, "X" : X, "F" : f_scf, "C" : c_scf, "P" : p_temp})
print("Электронная энергия:", ener_elec(p_temp, h_scf, f_scf))
print("Полная энергия:", ener_tot(r, Z, ener_elec(p_temp, h_scf, f_scf)))
return("Complete")
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}
