Python密度泛函理论

  • 1. 密度泛函理论
  • 2. Numpy 1D-DFT
    • 2.1 Differential operator
      • 2.1.1 First order differentiation
      • 2.1.2 Second oorder differentiation
    • 2.2 Non-interacting electrons
    • 2.3 Harmonic oscillator
    • 2.4 Density
    • 2.5 Exchange energy
    • 2.6 coulomb potential
    • 2.7 Solve the KS equation：Self-consistency loop

1. 密度泛函理论

密度泛函理论（Density Functional Theory）是一种研究多电子体系电子结构的量子力学方法。密度泛函理论在物理和化学上都有广泛的应用，特别是用来研究分子和凝聚态的性质，是凝聚态物理和计算化学领域最常用的方法之一。

电子结构理论的经典方法，特别是Hartree-Fock方法和后Hartree-Fock方法，是基于复杂的多电子波函数的。密度泛函理论的主要目标就是用电子密度取代波函数做为研究的基本量。因为多电子波函数有3N个变量（N 为电子数，每个电子包含三个空间变量），而电子密度仅是三个变量的函数，无论在概念上还是实际上都更方便处理。

Hohenberg-Kohn定理为泛函理论提供了坚实的理论基础：

  1. 第一定理指出体系的基态能量仅仅是电子密度的泛函。
  2. 第二定理证明了以基态密度为变量，将体系能量通过变分得到最小值之后就得到了基态能量。

密度泛函理论最普遍的应用是通过Kohn-Sham方法实现的。在Kohn-Sham DFT框架中，最难处理的多体问题被简化成了一个没有相互作用的电子在有效势场中的运动问题。

2. Numpy 1D-DFT

http://dcwww.camd.dtu.dk/~askhl/files/python-dft-exercises.pdf

Goal: write our own Kohn-Sham (KS) DFT code

Target: a harmonic oscillator including kinetic energy, electrostatic repulsion between the electrons, and the local density approximation for electronic interactions, ignoring correlation.

Hamiltonian:
H ^ = − 1 2 d 2 d x 2 + v ( x ) v ( x ) = v H a ( x ) + v L D A ( x ) + x 2 \hat{H}=-\frac{1}{2}\frac{d^2}{dx^2}+v(x)\\ v(x)=v_{Ha}(x)+v_{LDA}(x)+x^2 H^=−21​dx2d2​+v(x)v(x)=vHa​(x)+vLDA​(x)+x2

What we have to do?

  1. Represent the Hamiltonian
  2. Calculate the KS wavefunctions, the density
import numpy as np
import matplotlib.pyplot as plt

2.1 Differential operator

In order to represent kinetic operator

2.1.1 First order differentiation

approximate:

( d y d x ) i = y i + 1 − y i h (\frac{dy}{dx})_{i}=\frac{y_{i+1}-{y_{i}}}{h} (dxdy​)i​=hyi+1​−yi​​

then:

D i j = δ i + 1 , j − δ i , j h D_{ij}=\frac{\delta_{i+1,j}-\delta_{i,j}}{h} Dij​=hδi+1,j​−δi,j​​

we could write as follows:

( d y d x ) i = D i j y j (\frac{dy}{dx})_{i}=D_{ij} y_{j} (dxdy​)i​=Dij​yj​

The derivative may not be well defined at the end of the grid.

δ i j \delta_{ij} δij​ is Kronecker delta

Einstein summation is used for last equation

n_grid=200
x=np.linspace(-5,5,n_grid)
h=x[1]-x[0]
D=-np.eye(n_grid)+np.diagflat(np.ones(n_grid-1),1)
D = D / h

2.1.2 Second oorder differentiation

In the same way as the first order:

D i j 2 = δ i + 1 , j − 2 δ i , j + δ i − 1 , j h 2 D^2_{ij}=\frac{\delta_{i+1,j}-2\delta_{i,j}+\delta_{i-1,j}}{h^2} Dij2​=h2δi+1,j​−2δi,j​+δi−1,j​​

This could be written with the first order D i j D_{ij} Dij​, as follows (take care of the transpose):

D i j 2 = − D i k D j k D^2_{ij}=-D_{ik}D_{jk} Dij2​=−Dik​Djk​

The derivative may not be well defined at the end of the grid.

D2=D.dot(-D.T)
D2[-1,-1]=D2[0,0]

2.2 Non-interacting electrons

This is the Hamiltonian of non-interacting free particles in a box given by the size of grid:

H ^ = T ^ = − 1 2 d 2 d x 2 \hat{H} = \hat{T} = - \frac{1}{2} \frac{d^2}{dx^2} H^=T^=−21​dx2d2​

We could solve the KS equation as follows:

eig_non, psi_non=np.linalg.eigh(-D2/2)

plot (energies are shown in the label)

for i in range(3):
    plt.plot(x,psi_non[:,i], label=f"{eig_non[i]:.4f}")
    plt.legend(loc=1)

2.3 Harmonic oscillator

include the external potential v e x t = x 2 v_{ext}=x^2 vext​=x2:
H ^ = T ^ = − 1 2 d 2 d x 2 + x 2 \hat{H} = \hat{T} = - \frac{1}{2} \frac{d^2}{dx^2} + x^2 H^=T^=−21​dx2d2​+x2

we can write the potential as a matrix X X X, as follows:

X=np.diagflat(x*x)

and solve the KS.

eig_harm, psi_harm = np.linalg.eigh(-D2/2+X)

plot

for i in range(5):
    plt.plot(x,psi_harm[:,i], label=f"{eig_harm[i]:.4f}")
    plt.legend(loc=1)

2.4 Density

We will want to include the Coulomb or Hatree interacion as well as LDA exchange

Both of which are density functinals

So we need to calculate the electron density

Each state should be normalized:
∫ ∣ ψ ∣ 2 d x = 1 \int \lvert \psi \rvert ^2 dx = 1 ∫∣ψ∣2dx=1

let f n f_n fn​ be occupation numbers, the density n ( x ) n(x) n(x) can be written as follows:
n ( x ) = ∑ n f n ∣ ψ ( x ) ∣ 2 n(x)=\sum_n f_n \lvert \psi(x) \rvert ^2 n(x)=n∑​fn​∣ψ(x)∣2

Note:

  1. Each state fits up to two electrons: one with spin up, and one with spin down.
  2. In DFT, we calculate the ground state.
# integral
def integral(x,y,axis=0):
    dx=x[1]-x[0]
    return np.sum(y*dx, axis=axis)

number of electrons

num_electron=17

density

def get_nx(num_electron, psi, x):
    # normalization
    I=integral(x,psi**2,axis=0)
    normed_psi=psi/np.sqrt(I)[None, :]
    
    # occupation num
    fn=[2 for _ in range(num_electron//2)]
    if num_electron % 2:
        fn.append(1)

    # density
    res=np.zeros_like(normed_psi[:,0])
    for ne, psi  in zip(fn,normed_psi.T):
        res += ne*(psi**2)
    return res

plot

plt.plot(get_nx(num_electron,psi_non, x), label="non")
plt.plot(get_nx(num_electron,psi_harm, x), label="harm")
plt.legend(loc=1)

2.5 Exchange energy

Consider the exchange functional in the LDA and ignore the correlation for simplicity:

E X L D A [ n ] = − 3 4 ( 3 π ) 1 / 3 ∫ n 4 / 3 d x E_X^{LDA}[n] = -\frac{3}{4} \left(\frac{3}{\pi}\right)^{1/3} \int n^{4/3} dx EXLDA​[n]=−43​(π3​)1/3∫n4/3dx

The potential is given by the derivative of the exchange energy with respect to the density:

v X L D A [ n ] = ∂ E X L D A ∂ n = − ( 3 π ) 1 / 3 n 1 / 3 v_X^{LDA}[n] = \frac{\partial E_X^{LDA}}{\partial n} = - \left(\frac{3}{\pi}\right)^{1/3} n^{1/3} vXLDA​[n]=∂n∂EXLDA​​=−(π3​)1/3n1/3

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

2.6 coulomb potential

Electrostatic energy or Hatree energy

The expression of 3D-Hatree energy is not converged in 1D.

Hence we cheat and use a modified as follows:

E H a = 1 2 ∬ n ( x ) n ( x ′ ) ( x − x ′ ) 2 + ε d x d x ′ E_{Ha}=\frac{1}{2}\iint \frac{n(x)n(x')}{\sqrt{(x-x')^2+\varepsilon}}dxdx' EHa​=21​∬(x−x′)2+ε ​n(x)n(x′)​dxdx′

where ε \varepsilon ε is a small positive constant

The potential is given by:
v H a = ∫ n ( x ′ ) ( x − x ′ ) 2 + ε d x ′ v_{Ha}=\int \frac{n(x')}{\sqrt{(x-x')^2+\varepsilon}}dx' vHa​=∫(x−x′)2+ε ​n(x′)​dx′

In a matirx expression:
E H a = 1 2 n i n j h 2 ( x i − x j ) 2 + ε E_{Ha} = \frac{1}{2} \frac{n_in_jh^2}{\sqrt{(x_{i}-x_{j})^2+\varepsilon}} EHa​=21​(xi​−xj​)2+ε ​ni​nj​h2​
v H a , i = n j h ( x i − x j ) 2 + ε v_{Ha, i} = \frac{n_jh}{\sqrt{(x_{i}-x_{j})^2+\varepsilon}} vHa,i​=(xi​−xj​)2+ε ​nj​h​

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

2.7 Solve the KS equation：Self-consistency loop

  1. initialize the density (you can take an arbitrary constant)
  2. Calculate the Exchange and Hatree potentials
  3. Calculate the Hamiltonian
  4. Calculate the wavefunctions and eigen values
  5. If not converged, calculate the density and back to 1.
def print_log(i,log):
    print(f"step: {i:<5} energy: {log['energy'][-1]:<10.4f} energy_diff: {log['energy_diff'][-1]:.10f}")

max_iter=1000
energy_tolerance=1e-5
log={"energy":[float("inf")], "energy_diff":[float("inf")]}

nx=np.zeros(n_grid)
for i in range(max_iter):
    ex_energy, ex_potential=get_exchange(nx,x)
    ha_energy, ha_potential=get_hatree(nx,x)
    
    # Hamiltonian
    H=-D2/2+np.diagflat(ex_potential+ha_potential+x*x)
    
    energy, psi= np.linalg.eigh(H)
    
    # log
    log["energy"].append(energy[0])
    energy_diff=energy[0]-log["energy"][-2]
    log["energy_diff"].append(energy_diff)
    print_log(i,log)
    
    # convergence
    if abs(energy_diff) < energy_tolerance:
        print("converged!")
        break
    
    # update density
    nx=get_nx(num_electron,psi,x)
else:
    print("not converged")

plot

for i in range(5):
    plt.plot(x,psi[:,i], label=f"{energy[i]:.4f}")
    plt.legend(loc=1)

compare the density to free particles

plt.plot(nx)
plt.plot(get_nx(num_electron,psi_harm,x), label="no-interaction")
plt.legend()

参考文献来自桑鸿乾老师的课件：科学计算和人工智能
特此鸣谢！！！

更多推荐

【Python密度泛函理论】