用python实现密度泛函理论（1维）

gaohaikuo1

用python实现密度泛函理论https://blog.csdn.net/vor234/article/details/125088386

求解的总方程如下
$$\hat{\mathrm{H}}=-\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{dx}^2}+\mathrm{v}_{\mathrm{Ha}}\left(\mathrm{x}\right)+\mathrm{v}_{\mathrm{LDA}}\left(\mathrm{x}\right)+\mathrm{x}^2 $$
(1)不考虑势能时
$$\hat{\mathrm{H}}=\hat{\mathrm{T}}=-\frac{1}{2}\frac{\mathrm{d}^{2}}{\mathrm{dx}^{2}}$$
求解的关键在于如何处理二次微分,如果定义$(\frac{dy}{dx})_i=\frac{y_{i+1}-y_{i}}{h} $和$\mathrm{D}_{\mathrm{ij}}={\frac{\delta_{i+1,j}-\delta_{i,j}}{\mathrm{h}}} $（这里的i对应着矩阵的行，表示第几个点的微分；j表示列，表示每个点微分的计算过程），那么可以获得如下结论：
$$(\frac{\mathrm{dy}}{\mathrm{dx}})_i=\mathrm{D}_{ij}\mathrm{y}_j $$
这时，这个矩阵和y列向量相乘，就相当于后面的点减去了前面的点。
如果求两次导，则形式变为
$$\mathrm D_{ij}^2=\frac{\delta_{i+1,j}-2\delta_{i,j}+\delta_{i-1,j}}{\mathrm h^2}$$
最终的结果是两边之和减去2倍中间。
这时的代码如下

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. n_grid = 200
   
4. x = np.linspace(-5,5,n_grid)
   
5. h = x[1]-x[0]
   
6. D = -np.eye(n_grid)+np.diagflat(np.ones(n_grid-1),1)
   
7. D = D / h
   
8. D2=D.dot(-D.T)
   
9. D2[-1,-1]=D2[0,0]
   
10. eig_non, psi_non=np.linalg.eigh(-D2/2)
    
11. for i in range(3):
    
12.     print(i)   
    
13.     print(eig_non[i])
    
14.     plt.plot(x,psi_non[:,i], label=f"{eig_non[i]:.4f}")
    
15.     plt.legend(loc=1)
    
16.     plt.pause(0.01)
    

复制代码

gaohaikuo1

（2）考虑外部势能项$x^2$
需要注意的是，把该项要放在矩阵的对角线上，这是因为H矩阵的对角线是和当前点相关的，而其他地方的值指的是点之间的交互。

1. X=np.diagflat(x*x)
   
2. eig_harm, psi_harm = np.linalg.eigh(-D2/2+X)
   
3. for i in range(5):
   
4.     plt.plot(x,psi_harm[:,i], label=f"{eig_harm[i]:.4f}")
   
5.     plt.legend(loc=1)
   
6.     plt.pause(0.01)

复制代码

gaohaikuo1

（3）电子密度计算
接下来我们要计算Hatree和项和LDA交换项，然而这两个需要电子密度，因此需要先建立电子密度
电子密度为波函数的平方和

1. # integral
   
2. def integral(x,y,axis=0):
   
3.     dx=x[1]-x[0]
   
4.     return np.sum(y*dx, axis=axis)
   
5. 
   
6. num_electron=17
   
7. 
   
8. def get_nx(num_electron, psi, x):
   
9.     # normalization
   
10.     I=integral(x,psi**2,axis=0)
    
11.     normed_psi=psi/np.sqrt(I)[None, :]
    
12.    
    
13.     # occupation num
    
14.     fn=[2 for _ in range(num_electron//2)]
    
15.     if num_electron % 2:
    
16.         fn.append(1)
    
17. 
    
18.     # density
    
19.     res=np.zeros_like(normed_psi[:,0])
    
20.     for ne, psi  in zip(fn,normed_psi.T):
    
21.         res += ne*(psi**2)
    
22.     return res
    
23. 
    
24. plt.plot(get_nx(num_electron,psi_harm, x), label="harm")
    
25. plt.legend(loc=1)
    
26. plt.pause(0.01)

复制代码

gaohaikuo1

（3）求ex项LDA
根据密度可以求得LDA（忽略校正项）
$${E}_{X}^{LDA}[\mathrm{n}]=-\frac{3}{4}(\frac{3}{\pi})^{1/3}\int{n}^{4/3}\mathrm{dx}$$
$$\mathrm{v}_X^{\mathrm{LDA}}[\mathrm{n}]=\frac{\partial\mathrm{E}_X^{\mathrm{LDA}}}{\partial\mathrm{n}}=-(\frac{3}{\pi})^{1/3}\mathrm{n}^{1/3}$$

1. def get_exchange(nx,x):
   
2.     energy=-3./4.*(3./np.pi)**(1./3.)*integral(x,nx**(4./3.))
   
3.     potential=-(3./np.pi)**(1./3.)*nx**(1./3.)
   
4.     return energy, potential
   

复制代码

gaohaikuo1

(4)计算coulomb势
3D的Hatree在1D不收敛，因此用了一个近似。

1. def get_hatree(nx,x, eps=1e-1):
   
2.     h=x[1]-x[0]
   
3.     energy=np.sum(nx[None,:]*nx[:,None]*h**2/np.sqrt((x[None,:]-x[:,None])**2+eps)/2)
   
4.     potential=np.sum(nx[None,:]*h/np.sqrt((x[None,:]-x[:,None])**2+eps),axis=-1)
   
5.     return energy, potential

复制代码

这个实际上用了Numpy很有意思的一个算法
x[None,:]和x[:,None]实际是给x加了两个维度，比如原来x是10个元素的一维向量，x[None,:]就变为了1*10，而x[:,None]就变成了10*1，分别代表了xj和xi
换句话说，加入空维度None后，xi和xj在位置上就可以有区分，可以理解为冒号的位置，xj在第二个位置，xi在第一个位置
从公式也可以看到，最终的积分是对最后一个位置进行积分（axis = -1)，也就是对j进行积分。

gaohaikuo1

最终代码

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. n_grid = 200
   
4. x = np.linspace(-5,5,n_grid)
   
5. h = x[1]-x[0]
   
6. D = -np.eye(n_grid)+np.diagflat(np.ones(n_grid-1),1)
   
7. D = D / h
   
8. D2=D.dot(-D.T)
   
9. D2[-1,-1]=D2[0,0]
   
10. 
    
11. X=np.diagflat(x*x)
    
12. eig_harm, psi_harm = np.linalg.eigh(-D2/2+X)
    
13. 
    
14. 
    
15. # integral
    
16. def integral(x,y,axis=0):
    
17.     dx=x[1]-x[0]
    
18.     return np.sum(y*dx, axis=axis)
    
19. 
    
20. num_electron=17
    
21. 
    
22. def get_nx(num_electron, psi, x):
    
23.     # normalization
    
24.     I=integral(x,psi**2,axis=0)
    
25.     normed_psi=psi/np.sqrt(I)[None, :]
    
26.    
    
27.     # occupation num
    
28.     fn=[2 for _ in range(num_electron//2)]
    
29.     if num_electron % 2:
    
30.         fn.append(1)
    
31. 
    
32.     # density
    
33.     res=np.zeros_like(normed_psi[:,0])
    
34.     for ne, psi  in zip(fn,normed_psi.T):
    
35.         res += ne*(psi**2)
    
36.     return res
    
37. 
    
38. def get_exchange(nx,x):
    
39.     energy=-3./4.*(3./np.pi)**(1./3.)*integral(x,nx**(4./3.))
    
40.     potential=-(3./np.pi)**(1./3.)*nx**(1./3.)
    
41.     return energy, potential
    
42. 
    
43. def get_hatree(nx,x, eps=1e-1):
    
44.     h=x[1]-x[0]
    
45.     energy=np.sum(nx[None,:]*nx[:,None]*h**2/np.sqrt((x[None,:]-x[:,None])**2+eps)/2)
    
46.     potential=np.sum(nx[None,:]*h/np.sqrt((x[None,:]-x[:,None])**2+eps),axis=-1)
    
47.     return energy, potential
    
48. 
    
49. def print_log(i,log):
    
50.     print(f"step: {i:<5} energy: {log['energy'][-1]:<10.4f} energy_diff: {log['energy_diff'][-1]:.10f}")
    
51. 
    
52. max_iter=1000
    
53. energy_tolerance=1e-5
    
54. log={"energy":[float("inf")], "energy_diff":[float("inf")]}
    
55. 
    
56. nx=np.zeros(n_grid)
    
57. for i in range(max_iter):
    
58.     ex_energy, ex_potential=get_exchange(nx,x)
    
59.     ha_energy, ha_potential=get_hatree(nx,x)
    
60.    
    
61.     # Hamiltonian
    
62.     H=-D2/2+np.diagflat(ex_potential+ha_potential+x*x)
    
63.    
    
64.     energy, psi= np.linalg.eigh(H)
    
65.    
    
66.     # log
    
67.     log["energy"].append(energy[0])
    
68.     energy_diff=energy[0]-log["energy"][-2]
    
69.     log["energy_diff"].append(energy_diff)
    
70.     print_log(i,log)
    
71.    
    
72.     # convergence
    
73.     if abs(energy_diff) < energy_tolerance:
    
74.         print("converged!")
    
75.         break
    
76.    
    
77.     # update density
    
78.     nx=get_nx(num_electron,psi,x)
    
79. else:
    
80.     print("not converged")
    
81. 
    
82. plt.plot(nx)
    
83. plt.plot(get_nx(num_electron,psi_harm,x), label="no-interaction")
    
84. plt.legend()
    
85. plt.pause(0.01)
    

复制代码