zoukankan      html  css  js  c++  java
  • python实现: VMC做一维谐振子

    参考文献:J. M. Thijssen, 《Computational Physics》

    1. Formalism

    为了简化,设一维谐振子的哈密顿量为

    \[\hat{H} = - \frac{d^2}{dx^2} + x^2 \]

    引入拟设波函数

    \[\psi(x, \alpha) = e^{-\alpha x^2}, \]

    定义局域能量

    \[E_L(x, \alpha) = \frac{ \hat{H} \psi(x, \alpha) }{ \psi(x, \alpha)}, \]

    如果 \(\alpha\) 的值恰当,使得 \(\psi(x,\alpha)\)\(\hat{H}\) 的本征态,则 \(\forall x, E_L(x,\alpha) = \epsilon\)
    变分法中的“平均能量”为

    \[\langle E(\alpha) \rangle = \int E_L(x, \alpha) \psi(x,\alpha)^2 dx. \]

    所以,设置大量的随机行者,用 Markov-Metropolis 方法,使得行者的分布密度 \(\rho(x,\alpha) \propto \psi(x,\alpha)^2\),然后做 \(E_L(x,\alpha)\)的抽样平均,即得 \(\langle E(\alpha) \rangle\) 的近似值。
    另外,可以观察抽样方差 \(\sigma^2(\alpha) = \langle E(\alpha)^2 \rangle - \langle E(\alpha) \rangle^2\)\(\alpha\) 的函数关系。如果 \(\alpha\) 恰当,使得 \(\psi(x,\alpha)\)\(\hat{H}\) 的本征态,则 \(\sigma^2(\alpha) = 0\),如果 \(\alpha\) 不够理想,就有 \(\sigma^2(\alpha) > 0\)

    2. 代码

    
    # 一维谐振子求基态能量
    import numpy as np
    import matplotlib.pyplot as plt
    
    # trial wave: psi = e^{ -alpha x^2 }
    # H = - d^2/dx^2 + x*x
    import scipy.integrate
    
    def psi(x,alpha):
        return np.exp( -alpha * x * x )
    
    #help(scipy.integrate.quad); exit(1)
    def plotpsi2(alpha):
        C = scipy.integrate.quad(psi,-np.inf,np.inf,args=(2*alpha,))[0]
        x = np.arange(-3,3,0.1)
        y = [  psi(t,alpha) * psi(t,alpha) / C for t in x ]
        plt.plot(x,y)
    #plotpsi2(0.4); exit(1)
    
    # d/dx psi = -2 alpha x e^{- alpha * x * x }
    # d^2 /d x^2 psi = -2 alpha + 4 * alpha^2 x^2 e^{-alpha x^2}
    # ( - d^2/dx^2 + x*x ) psi / psi =  2 alpha + (- 4 alpha^2 + 1) x^2
    def EL( alpha, x):
        return 2*alpha + (1-4*alpha*alpha)*x*x
    
    # if alpha = 0.5, EL(alpha, x) = 1
    # if alpha = 0.6, EL(alpha, x) = 1.2 - 0.44 * x*x
    
    def Accept( x1, x2, alpha, psi ):
        if( abs(x2)>3 ): return 0
        p = ( psi(x2,alpha) / psi(x1,alpha) )**2
        #print("p=",p)
        if p>1: return 1
        else: return p
    
    def tryonestep(x,nwalker,h,alpha):
        for j in range(nwalker):
            step = h * ( 2*np.random.randint(0,2) - 1 )
            A = Accept(x[j], x[j] + step, alpha, psi)
            if (np.random.random() < A):
                x[j] += step
    
    def ELsampling( alpha, nwalker, psi, EL, a, b, h, nstep ):
        x = (b-a)*np.random.random(nwalker) + a # 初始位置
        for i in range(nstep):
            tryonestep(x, nwalker, h, alpha)
        #plotpsi2(alpha);
        #plt.hist(x, bins=np.arange(-3.1,3.1,0.2), rwidth=0.8, density="True" ); plt.show()
        #print( "alpha = ",alpha, " <E>_L = ", np.average( [ EL(alpha,t) for t in x ] ) )
        return [ EL(alpha,t) for t in x ]
    
    aveE = []; aveE2 = [] # < E >, < E^2 >
    Alpha = np.arange(0.4,0.6,0.02)
    for alpha in Alpha:
        ELsample = ELsampling(alpha, 1000, psi, EL, a=-3, b=3, h=0.2, nstep=100 )
        aveE.append( np.average(ELsample) )
        aveE2.append( np.average([t*t for t in ELsample ]) )
    sigma2 = aveE2 - np.array([t*t for t in aveE])
    plt.plot(Alpha, sigma2)
    plt.show()
    

    3. 结果

    做出来 \(\sigma^2(\alpha) \sim \alpha\) 的图如下。可见 \(\alpha=0.5\) 是最优参数。
    image

    4. 小结

    • 代码中有绘制 \(\rho(x,\alpha)\) 的片段,可以用于观察随机行者在接受概率的指引下收敛于 \(\rho(x,\alpha)\) 的速度。在我的试验中,行者处在 \(-3, -2.8, \cdots, 2.8, 3.0\) 这些格点上,一共大约 30 个位置,大约 100 步随机行走后,行者分布就已经收敛至 \(\rho(x, \alpha)\) 曲线。
    • \(\langle E(\alpha) \rangle\) 远远没有 \(\sigma^2(\alpha)\) 犀利。上面的代码以及其中的参数跑出来的 \(\langle E(\alpha) \rangle\) 曲线非常磕碜,但是上面展示的 \(\sigma^2(\alpha)\) 曲线就非常规整,且容易反映出问题。如果你开个根号,画 \(\sigma \sim \alpha\),就会更加 sharp。
  • 相关阅读:
    85. Maximal Rectangle
    120. Triangle
    72. Edit Distance
    39. Combination Sum
    44. Wildcard Matching
    138. Copy List with Random Pointer
    91. Decode Ways
    142. Linked List Cycle II
    异或的性质及应用
    64. Minimum Path Sum
  • 原文地址:https://www.cnblogs.com/luyi07/p/15621162.html
Copyright © 2011-2022 走看看