import hfscf
import numpy as np
%pip install sympy
/usr/lib/python3/dist-packages/secretstorage/dhcrypto.py:15: CryptographyDeprecationWarning: int_from_bytes is deprecated, use int.from_bytes instead
from cryptography.utils import int_from_bytes
/usr/lib/python3/dist-packages/secretstorage/util.py:19: CryptographyDeprecationWarning: int_from_bytes is deprecated, use int.from_bytes instead
from cryptography.utils import int_from_bytes
Defaulting to user installation because normal site-packages is not writeable
Collecting sympy
Downloading sympy-1.7.1-py3-none-any.whl (5.9 MB)
|████████████████████████████████| 5.9 MB 2.9 MB/s eta 0:00:01
Collecting mpmath>=0.19
Downloading mpmath-1.2.1-py3-none-any.whl (532 kB)
|████████████████████████████████| 532 kB 33.6 MB/s eta 0:00:01
Installing collected packages: mpmath, sympy
WARNING: The script isympy is installed in '/home/students/y17/skorpion7/.local/bin' which is not on PATH.
Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
Successfully installed mpmath-1.2.1 sympy-1.7.1
Note: you may need to restart the kernel to use updated packages.
def SCF (r = 1.4632, Z=[1,1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3, vbs=False):
# Нужно задать центры атомов
R = [0, r]
# Нужно сгенерировать матрицу перекрывания (S), de traslape
s_scf = hfscf.S(R, b1, b2, b)
# Нужно сгенерировать гамильтониан (H)
h_scf = hfscf.H(R, Z, b1, b2, b)
# Нужно диагонализировать матрицу S и найти матрицу X. Здесь Xa - это сопряженно-транспонированная матрица
X = hfscf.diagon(m=s_scf)
Xa = X.getH()
# Нужно чоздать матрицу плотности P
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64)
# Начинается проверка сходимсоти F. Итеративно заполняется матрица F, диагонализируется, проверяется разница
for iteracion in range(50):
# Построение матрицы Фока как F = H + G. Перед этим - расчёт двух электронных частей матрицы Фока (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_scf = X * c_tra
# Теперь из матрицы С нужно построить матрицу P. Пересчитать матрицу плотности P
p_temp = hfscf.P(C=c_scf)
print("Завершение " + str(iteracion + 1) + "итерации")
# Теперь проверка сходимости
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}
# Проверим запуск нашей функции в режиме параметров по умолчанию
prowerochka = SCF()
Завершение 1итерации Завершение 2итерации Завершение 3итерации Завершение 4итерации Завершение 5итерации Завершение 6итерации Завершение 7итерации SCF сошлась!
# На выходе получаются посчитанные матрицы
prowerochka
{'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 ]])}
# Теперь можно нарисовать орбитали в первом измерении
hfscf.orbital(prowerochka['C'])
# Также их можно нарисовать и во втором измерении
hfscf.orbital2D(prowerochka['C'], prowerochka['F'], prowerochka['X'])
# То что слева - связывающая орбиталь, то что слева - разрыхляющая орбиталь
# Теперь расчёт энергий
Energías_orbitales = hfscf.ener_orbs(prowerochka['X'], prowerochka['F'])
Energía_electrónica_de_un_molécula = hfscf.ener_elec(prowerochka['P'], prowerochka['H'], prowerochka['F'])
energía_total = hfscf.ener_tot()
print('Энергии орбиталей: ' + str(Energías_orbitales[0]) + ' и ' + str(Energías_orbitales[1]))
Энергии орбиталей: -0.5646153594583119 и 0.6368673980935621
print('Электронная энергия молекулы: ' + str(Energía_electrónica_de_un_molécula))
Электронная энергия молекулы: -1.7974485548087504
print('Полная энергия молекулы: ' + str(energía_total))
Полная энергия молекулы: 0.683433570256971