import numpy as np
from scipy.linalg import eigh
from math import pi
import matplotlib.pyplot as plt
a = [13.00773, 1.962079, 0.444529, 0.1219492]
S = np.zeros([4,4])
T = np.zeros([4,4])
A = np.zeros([4,4])
H = np.zeros([4,4])
for i in range(0, 4):
for j in range(0,4):
S[i][j] = (pi /(a[i]+a[j]))**(3/2)
T[i][j] = 3 * a[i]* a[j] * pi**(3/2) /(a[i]+ a[j])**(5/2)
A[i][j] = -2*pi / (a[i] + a[j])
H = T + A
eigvals, eigvecs = eigh(H, S)
print(eigvals.min())
参考计算物理J.M.Thijssen P36页内容，涉及到一个广义特征值的求解
