import numpy as np
import math
import scipy.special
import scipy.misc
import hfscf
def SCF (r = 1.4632, Z=[1,1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3, vbs=False):
#Задаем координты ядер
R = [0, r]
#Создаем матрицу перекрывания S
s_scf = hfscf.S(R, b1, b2, b)
#Создаем гамильтониан Н
h_scf = hfscf.H(R, Z, b1, b2, b)
#Диагонализуем S, получая Х
X = hfscf.diagon(s_scf)
Xa = X.getH()
#Получаем матрицу рассчитанной плотности Р
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64) # Referencia (7) p. 148
#Начинаем итерационное улучшение
for iteracion in range(50):
#Получаем матрицу Фока
# F = H + G
g_scf = hfscf.G(r, p_scf, b1, b2, b)
f_scf = h_scf + g_scf # Referencia (7) p. 141 eq. (3.154)
#Создаем матрицу F'
f_tra = Xa * f_scf * X
#Диагонализуем матрицу F, получая матрицу C'
c_tra = hfscf.diagon2(f_tra)
#Получаем матрицу коэффициентов C
c_scf = X * c_tra
#Получаем матрицу плотности P на основе C
p_temp = hfscf.P(c_scf)
#Выводим энергии
electron_E = hfscf.ener_elec(p_temp, h_scf, f_scf)
total_E = hfscf.ener_tot(r, Z, 0)
print("Энергия электронов: ", electron_E, ', общая энергия: ', total_E)
#Рисуем промежуточные картиночки
# hfscf.orbital(c_tra, r, b1, b2, b)
# hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, delta=0.02)
print("Конец " + str(iteracion + 1) + " итерации.")
print('='*70)
#Проверяем сходимость и делаем выводы
if np.linalg.norm(p_temp - p_scf) < 1E-4:
print('='*70)
print('='*70)
print("\n\n-->Поле сошлось!")
print("Итоговая энергия электронов: \t", electron_E, '\nИтоговая общая энергия: \t', total_E)
orb_E = hfscf.ener_orbs(X, f_scf)
print('Энергии орбиталей: ' + str(orb_E[0]) + ' и ' + str(orb_E[1]))
hfscf.orbital(c_tra, r, b1, b2, b)
hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, delta=0.02)
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()
Энергия электронов: -1.72842063181647 , общая энергия: 0.683433570256971 Конец 1 итерации. ====================================================================== Энергия электронов: -1.565429982099157 , общая энергия: 0.683433570256971 Конец 2 итерации. ====================================================================== Энергия электронов: -1.795075334555484 , общая энергия: 0.683433570256971 Конец 3 итерации. ====================================================================== Энергия электронов: -1.7974292055137373 , общая энергия: 0.683433570256971 Конец 4 итерации. ====================================================================== Энергия электронов: -1.7974483974402236 , общая энергия: 0.683433570256971 Конец 5 итерации. ====================================================================== Энергия электронов: -1.7974485535391342 , общая энергия: 0.683433570256971 Конец 6 итерации. ====================================================================== Энергия электронов: -1.7974485548087507 , общая энергия: 0.683433570256971 Конец 7 итерации. ====================================================================== ====================================================================== ====================================================================== -->Поле сошлось! Итоговая энергия электронов: -1.7974485548087507 Итоговая общая энергия: 0.683433570256971 Энергии орбиталей: -0.5646153594583119 и 0.6368673980935621
{'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 ]])}