高斯函数展开计算氦原子基态能量

导入模块

import numpy as np
from scipy.linalg import eigh
from math import pi,sqrt

高斯基函数和动能，势能矩阵，哈密顿量矩阵

a = [0.298073, 1.242567, 5.782948, 38.474970]
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 + 2*A

电子电子库伦作用矩阵

Q = np.zeros([4,4,4,4])
for i in range(0,4):
    for j in range(0,4):
        for k in range(0,4):
            for l in range(0,4):
                Q[i][j][k][l] = 2*pi**(5/2)/((a[i]+a[j])*(a[k]+a[l])*sqrt(a[i]+a[j]+a[k]+a[l]))

初值波函数和F矩阵

C1 = np.ones([4,1])
F = np.zeros([4,4])

迭代过程

for i in range(0,100):
    C1 = C1.reshape(4,1)
    C2 = np.transpose(C1)
    c = np.dot(np.dot(C2,S),C1)
    C1 = C1/ sqrt(c)
    C2 = C2/ sqrt(c)
    C1 = C1.reshape(4,)
    temp = np.einsum('pqrs, s, r ->pq', Q, C1, C1)
    F = H + temp
    eigvals, eigvecs = eigh(F, S)
    C1 = eigvecs[:,0]

输出基态能量

C1 = C1.reshape(4,1)
C2 = np.transpose(C1)
c = np.dot(np.dot(C2,S),C1)
C1 = C1/ sqrt(c)
C2 = C2/ sqrt(c)
E1 = 2 * np.dot(C2,np.dot(H,C1))
C1 = C1.reshape(4,)
E2 =  np.einsum('pqrs, s, r,q,p ->', Q, C1, C1,C1,C1)
E = E1 + E2

参考计算物理J.M.Thijssen P51页内容，其中涉及的张量运算使用numpy.einsum函数