一个简单而不平凡的体系
一维双势阱模型
在 7 月 12 日的文章 中，我们讨论了非平衡分子动力学中的 Jarzynski 等式。我们先简单回顾一下 Jarzynski 等式适用的条件：
考虑一个经典、恒温的系统，它具有 N N N x = ( q , p ) x=(q,p) x = ( q , p )
现设系统经历了一个过程并改变了它的 Hamilton 量，则在这个过程中 Hamilton 量可以写作 H ( x , λ ) H(x,\lambda) H ( x , λ ) λ = 0 \lambda=0 λ = 0 λ = 1 \lambda=1 λ = 1 λ \lambda λ λ ( t ) \lambda(t) λ ( t )
取定初始时刻相空间坐标为 x ( 0 ) x(0) x ( 0 ) x ( t ) x(t) x ( t )
W [ x ( t ) ] = ∫ 0 τ ∂ H ( x , λ ) ∂ λ d λ ( t ) d t d t W[\mathbf x(t)]=\int_{0}^{\tau} \frac{\partial H(\mathbf x,\lambda)}{\partial \lambda} \frac{\mathrm d\lambda(t) }{\mathrm d t} \mathrm d t W [ x ( t )] = ∫ 0 τ ∂ λ ∂ H ( x , λ ) d t d λ ( t ) d t 现在我们给出一个这样的 Hamilton 量（尽管体系非常简单，它用于几乎所有文献的最初的测试）：
H λ = 1 2 p 2 + q 4 + 16 ( λ − 1 ) q 2 H_\lambda=\frac12p^2+q^4+16(\lambda-1)q^2 H λ = 2 1 p 2 + q 4 + 16 ( λ − 1 ) q 2 除去普通的动能项，剩下的势能项由 λ = 0 \lambda=0 λ = 0 q 4 − 16 q 2 q^4-16q^2 q 4 − 16 q 2 λ = 1 \lambda=1 λ = 1 q 4 q^4 q 4
给出具体的 Hamilton 量，我们只需要在恒温条件下求解 Hamilton 正则方程并计算相应的功即可。
含时 Hamilton 量的二阶辛算法
扩展 Hamilton 量和正则方程的构建
在 7 月 19 日的文章 中，我们总结了在 Hamilton 量不含时的情况下，辛算法的一般构建方式。不过上面的讨论要求我们处理一个 Hamilton 量随标度参量变化，因而也就随时间变化的情况：H = H ( p , q , λ ( t ) ) H=H(p,q,\lambda (t)) H = H ( p , q , λ ( t ))
先考虑一般的情况，即 R 2 n \mathbb R^{2n} R 2 n
我们引入虚拟位置 q n + 1 = t q_{n+1}=t q n + 1 = t H ( p 1 , ⋯ , p n , q 1 , ⋯ , q n + 1 H(p_1,\cdots,p_n,q_1,\cdots,q_{n+1} H ( p 1 , ⋯ , p n , q 1 , ⋯ , q n + 1 p n + 1 = − H p_{n+1}=-H p n + 1 = − H K = H + p n + 1 K=H+p_{n+1} K = H + p n + 1
q ˙ i = ∂ K ∂ p i p ˙ i = − ∂ K ∂ q i \begin{aligned}
\dot q_i&=\frac{\partial K}{\partial p_i}\\
\dot p_i&=-\frac{\partial K}{\partial q_i}\\
\end{aligned} q ˙ i p ˙ i = ∂ p i ∂ K = − ∂ q i ∂ K 而新加入的变量满足的方程是
q ˙ n + 1 = ∂ K ∂ p n + 1 = 1 p ˙ n + 1 = − ∂ K ∂ q n + 1 \begin{aligned}
\dot q_{n+1}&=\frac{\partial K}{\partial p_{n+1}}=1\\
\dot p_{n+1}&=-\frac{\partial K}{\partial q_{n+1}}\\
\end{aligned} q ˙ n + 1 p ˙ n + 1 = ∂ p n + 1 ∂ K = 1 = − ∂ q n + 1 ∂ K 可见这是自洽的。虽然 q n + 1 q_{n+1} q n + 1 t t t K K K 2 n + 2 2n+2 2 n + 2 R 2 n + 2 R^{2n+2} R 2 n + 2
Nosé-Poincaré 控温的实现
考虑上含时的势能之后（我们将时间记作是 Q Q Q
H = s [ p 2 2 m s 2 + V ( q , Q ) + ln s + π 2 2 M − H 0 ] H=s\left[\frac{p^2}{2ms^2}+V(q,Q)+\ln s+\frac{\pi^2}{2M}-H_0\right] H = s [ 2 m s 2 p 2 + V ( q , Q ) + ln s + 2 M π 2 − H 0 ] 以及升维之后的 K = H + P K=H+P K = H + P
H 1 = s [ p 2 2 m s 2 + ln s − H 0 ] H 2 = s V ( q ) H 3 = s π 2 2 M \begin{aligned}
H_1&=s\left[\frac{p^2}{2ms^2}+ \ln s-H_{0}\right]\\
H_2&=sV(q)\\
H_3&=\frac{s\pi^2}{2M}
\end{aligned} H 1 H 2 H 3 = s [ 2 m s 2 p 2 + ln s − H 0 ] = s V ( q ) = 2 M s π 2 我们发现，如果将 V ( q , Q ) V(q,Q) V ( q , Q ) V ( q ) V(q) V ( q ) H 2 H_2 H 2 H 2 H_2 H 2 P P P H 1 H_1 H 1 H 1 H_1 H 1 p , P p,P p , P s s s K = K 1 + K 2 + K 3 K=K_1+K_2+K_3 K = K 1 + K 2 + K 3
K 1 = P + s [ p 2 2 s 2 + ln s − H 0 ] K 2 = s V ( q , Q ) K 3 = s π 2 2 M \begin{aligned}
K_1&=P+s\left[\frac{p^2}{2s^2}+\ln s-H_{0}\right]\\
K_2&=sV(q,Q)\\
K_3&=\frac{s\pi^2}{2M}
\end{aligned} K 1 K 2 K 3 = P + s [ 2 s 2 p 2 + ln s − H 0 ] = s V ( q , Q ) = 2 M s π 2 对应的辛差分格式仍然适用：x n + 1 = e Δ t H 3 / 2 e Δ t H 2 / 2 e Δ t H 1 e Δ t H 2 / 2 e Δ t H 3 / 2 x n x_{n+1}=e^{\Delta tH_3/2}e^{\Delta tH_2/2}e^{\Delta tH_1}e^{\Delta tH_2/2}e^{\Delta tH_3/2}x_n x n + 1 = e Δ t H 3 /2 e Δ t H 2 /2 e Δ t H 1 e Δ t H 2 /2 e Δ t H 3 /2 x n x 是 6 维向量。我们具体将这个格式写出来（高能预警 ）：
Evolve H 3 / 2 : s ∗ = s n ( 1 + π n 2 M Δ t 2 ) 2 π ∗ = π n / ( 1 + π n 2 M Δ t 2 ) Evolve H 2 / 2 : p ∗ = p n − s ∗ ∂ V ∂ q ∣ q n , Q n Δ t 2 P ∗ = P − s ∗ ∂ V ∂ Q ∣ q n , Q n Δ t 2 π ∗ ∗ = π ∗ − V ( q n , Q n ) Δ t 2 Evolve H 1 : q n + 1 = q n + p ∗ s ∗ Δ t Q n + 1 = Q n + Δ t π ∗ ∗ ∗ = π ∗ ∗ + [ 1 2 ( p ∗ s ∗ ) 2 − ln s ∗ + H 0 − 1 ] Δ t Evolve H 2 / 2 : p n + 1 = p ∗ − s ∗ ∂ V ∂ q ∣ q n + 1 , Q n + 1 Δ t 2 P n + 1 = P ∗ − s ∗ ∂ V ∂ Q ∣ q n + 1 , Q n + 1 Δ t 2 π ∗ ∗ ∗ ∗ = π ∗ ∗ ∗ − V ( q n + 1 , Q n + 1 ) Δ t 2 Evolve H 3 / 2 : s n + 1 = s ∗ ( 1 + π ∗ ∗ ∗ ∗ 2 M Δ t 2 ) 2 π n + 1 = π ∗ ∗ ∗ ∗ / ( 1 + π ∗ ∗ ∗ ∗ 2 M Δ t 2 ) \begin{aligned}
&\text{Evolve}~H_3/2:\\
s^* &=s^{n}\left(1+\frac{\pi^{n}}{2M} \frac{\Delta t}{2}\right)^{2} \\
\pi^{*} &=\pi^{n} /\left(1+\frac{\pi^{n}}{2M} \frac{\Delta t}{2}\right) \\
&\text{Evolve}~H_2/2:\\
p^* &=p^{n}-\left.s^* \frac{\partial V}{\partial q}\right|_{q_n,Q_n} \frac{\Delta t}{2} \\
P^*&=P-\left.s^*\frac{\partial V}{\partial Q}\right|_{q_n,Q_n} \frac{\Delta t}{2}\\
\pi^{**}&=\pi^*-V(q_n,Q_n)\frac{\Delta t}{2} \\
&\text{Evolve}~H_1:\\
q^{n+1} &=q^{n}+\frac{p^*}{s^*} \Delta t \\
Q^{n+1} &=Q^{n}+ \Delta t \\
\pi^{***}&= \pi^{**}+\left[\frac{1}{2}\left(\frac{p^*}{s^*}\right)^{2} -\ln s^*+H_{0}-1\right] \Delta t \\
&\text{Evolve}~H_2/2:\\
p^{n+1} &=p^*-\left.s^* \frac{\partial V}{\partial q}\right|_{q_{n+1},Q_{n+1}} \frac{\Delta t}{2}\\
P^{n+1}&=P^*-\left.s^*\frac{\partial V}{\partial Q}\right|_{q_{n+1},Q_{n+1}} \frac{\Delta t}{2}\\
\pi^{****}&=\pi^{***}-V(q_{n+1},Q_{n+1})\frac{\Delta t}{2} \\
&\text{Evolve}~H_3/2:\\
s^{n+1}&= s^*\left(1+\frac{\pi^{****}}{2M} \frac{\Delta t}{2}\right)^{2} \\ \pi^{n+1}&= \pi^{*** *} /\left(1+\frac{\pi^{****}}{2M} \frac{\Delta t}{2}\right) \end{aligned} s ∗ π ∗ p ∗ P ∗ π ∗∗ q n + 1 Q n + 1 π ∗∗∗ p n + 1 P n + 1 π ∗∗∗∗ s n + 1 π n + 1 Evolve H 3 /2 : = s n ( 1 + 2 M π n 2 Δ t ) 2 = π n / ( 1 + 2 M π n 2 Δ t ) Evolve H 2 /2 : = p n − s ∗ ∂ q ∂ V q n , Q n 2 Δ t = P − s ∗ ∂ Q ∂ V q n , Q n 2 Δ t = π ∗ − V ( q n , Q n ) 2 Δ t Evolve H 1 : = q n + s ∗ p ∗ Δ t = Q n + Δ t = π ∗∗ + [ 2 1 ( s ∗ p ∗ ) 2 − ln s ∗ + H 0 − 1 ] Δ t Evolve H 2 /2 : = p ∗ − s ∗ ∂ q ∂ V q n + 1 , Q n + 1 2 Δ t = P ∗ − s ∗ ∂ Q ∂ V q n + 1 , Q n + 1 2 Δ t = π ∗∗∗ − V ( q n + 1 , Q n + 1 ) 2 Δ t Evolve H 3 /2 : = s ∗ ( 1 + 2 M π ∗∗∗∗ 2 Δ t ) 2 = π ∗∗∗∗ / ( 1 + 2 M π ∗∗∗∗ 2 Δ t ) 当然实际计算的过程中我们并不关心这个 Hamilton 量到底变化了多少，所以 P P P Q Q Q
误差理论
点估计理论
为了定量刻画不充分采样产生的误差，我们需要一点统计学的分析。为简便起见，移动能量零点使初态 F = 0 F=0 F = 0 F = ln ⟨ e − W ⟩ F=\ln\langle e^{-W}\rangle F = ln ⟨ e − W ⟩ N N N F N F_N F N
F N = − ln ( 1 N ∑ i = 1 N e − W i ) F_N=-\ln\left(\frac1N\sum_{i=1}^Ne^{- W_i}\right) F N = − ln ( N 1 i = 1 ∑ N e − W i ) 理想的估值方案需要在尽可能少的计算时间内获得尽可能精确的自由能，为此我们定义均方误差：ε 2 [ F N ] = ⟨ ( F N − F ) 2 ⟩ \varepsilon^2[F_N]=\langle(F_N-F)^2\rangle ε 2 [ F N ] = ⟨( F N − F ) 2 ⟩ ⟨ ( F N − F ) 2 ⟩ = ⟨ ( F N − ⟨ F N ⟩ ) 2 ⟩ + ( ⟨ F N ⟩ − F ) 2 \langle(F_N-F)^2\rangle=\langle(F_N-\langle F_N\rangle)^2\rangle+(\langle F_N\rangle-F)^2 ⟨( F N − F ) 2 ⟩ = ⟨( F N − ⟨ F N ⟩ ) 2 ⟩ + (⟨ F N ⟩ − F ) 2 σ 2 [ F N ] \sigma^2[F_N] σ 2 [ F N ] b 2 [ F N ] b^2[F_N] b 2 [ F N ] ε 2 [ F N ] = σ 2 [ F N ] + b 2 [ F N ] \varepsilon^2[F_N]=\sigma^2[F_N]+b^2[F_N] ε 2 [ F N ] = σ 2 [ F N ] + b 2 [ F N ]
直接指数平均是有偏估计
上述最简单的估值方案表现如何？首先考虑它的系统误差 b 2 b^2 b 2
X = e − W , X N = 1 N ∑ i = 1 N e − W i X=e^{-W},X_N=\frac1N\sum_{i=1}^Ne^{-W_i} X = e − W , X N = N 1 i = 1 ∑ N e − W i 也即 F N = − ln X N F_N=-\ln X_N F N = − ln X N ⟨ F N ⟩ = − ⟨ ln X N ⟩ \langle F_N\rangle=-\langle\ln X_N\rangle ⟨ F N ⟩ = − ⟨ ln X N ⟩
由 Jensen’s（琴生）不等式，ln ⟨ X N ⟩ < ⟨ ln X N ⟩ \ln\langle X_N\rangle<\langle\ln X_N\rangle ln ⟨ X N ⟩ < ⟨ ln X N ⟩ ⟨ X N ⟩ = ⟨ X ⟩ = e − F \langle X_N\rangle=\langle X\rangle=e^{-F} ⟨ X N ⟩ = ⟨ X ⟩ = e − F ⟨ F N ⟩ ≥ − ln e − F = F \langle F_N\rangle\ge -\ln e^{-F}=F ⟨ F N ⟩ ≥ − ln e − F = F
这说明我们总是高估真实的自由能变。
高估了多少？我们记随机变量 X X X m m m s 2 s^2 s 2 ⟨ X ⟩ = e − F = m \langle X\rangle=e^{-F}=m ⟨ X ⟩ = e − F = m ln ( ⋅ ) \ln(\cdot) ln ( ⋅ ) X N / m X_N/m X N / m
− ⟨ ln X N ⟩ = − ln m − ⟨ ln X N m ⟩ = − ln m − ⟨ Y N − Y N 2 2 + Y N 3 3 + . . . ⟩ -\langle \ln X_N\rangle=-\ln m-\left\langle\ln \frac{X_N}m\right\rangle=-\ln m-\left\langle Y_N-\frac{Y_N^2}2+\frac{Y_N^3}3+...\right\rangle − ⟨ ln X N ⟩ = − ln m − ⟨ ln m X N ⟩ = − ln m − ⟨ Y N − 2 Y N 2 + 3 Y N 3 + ... ⟩ 其中 Y N = X N / m − 1 Y_N=X_N/m-1 Y N = X N / m − 1 − ln m -\ln m − ln m ⟨ Y N ⟩ = 0 \langle Y_N\rangle=0 ⟨ Y N ⟩ = 0
b N ≈ 1 2 ⟨ Y N 2 ⟩ = 1 2 ⟨ ( X N − m ) 2 m 2 ⟩ = s 2 2 N m 2 b_N\approx \frac12\left\langle Y_N^2\right\rangle=\frac12\left\langle\frac{(X_N-m)^2}{m^2}\right\rangle=\frac{s^2}{2Nm^2} b N ≈ 2 1 ⟨ Y N 2 ⟩ = 2 1 ⟨ m 2 ( X N − m ) 2 ⟩ = 2 N m 2 s 2 同理，可以算出随机误差
σ N 2 ≈ ⟨ Y N 2 ⟩ = s 2 N m 2 \sigma_N^2\approx \langle Y_N^2\rangle=\frac{s^2}{Nm^2} σ N 2 ≈ ⟨ Y N 2 ⟩ = N m 2 s 2 鉴于 b N 2 = O ( N − 2 ) b_N^2=O(N^{-2}) b N 2 = O ( N − 2 ) σ N 2 \sigma_N^2 σ N 2 ε N 2 \varepsilon_N^2 ε N 2 ε N 2 ≈ σ N 2 \varepsilon_N^2\approx\sigma_N^2 ε N 2 ≈ σ N 2
又由于 X = e − W X=e^{-W} X = e − W ε N 2 \varepsilon_N^2 ε N 2 W d W_d W d W − F W-F W − F
ε N 2 ≈ ⟨ ( e − W d − 1 ) 2 ⟩ N \varepsilon_N^2\approx \frac{\langle (e^{-W_d}-1)^2\rangle}N ε N 2 ≈ N ⟨( e − W d − 1 ) 2 ⟩ 换言之，获得准确度为 k T kT k T N ≈ ⟨ ( e − W d − 1 ) 2 ⟩ N\approx\langle(e^{-W_d}−1)^2\rangle N ≈ ⟨( e − W d − 1 ) 2 ⟩ ⟨ ( e − W d − 1 ) 2 ⟩ \langle(e^{-W_d}−1)^2\rangle ⟨( e − W d − 1 ) 2 ⟩
一些结果
解析计算
自由能可以用配分函数表示：
F = − ln Z 1 Z 0 = − ln { ∫ d p d q exp [ − ( 1 2 p 2 + q 4 − 16 q 2 ) ] } = ln π e 32 ( I − 1 / 4 ( 32 ) + I 1 / 4 ( 32 ) ) 2 Γ ( 5 / 4 ) ≈ 62.9407 \begin{aligned}
F&=-\ln\frac{Z_1}{Z_0}\\
&=-\ln\left\{\int \mathrm dp\mathrm dq\exp\left[-\left(\frac12p^2+q^4-16q^2\right)\right]\right\}\\
&=\ln \frac{\pi e^{32}\left(I_{-1 / 4}(32)+I_{1 / 4}(32)\right)}{\sqrt{2} \Gamma(5 / 4)}\\
&\approx 62.9407
\end{aligned} F = − ln Z 0 Z 1 = − ln { ∫ d p d q exp [ − ( 2 1 p 2 + q 4 − 16 q 2 ) ] } = ln 2 Γ ( 5/4 ) π e 32 ( I − 1/4 ( 32 ) + I 1/4 ( 32 ) ) ≈ 62.9407 如果初态取自 Boltzmann 分布，容易得到当 λ \lambda λ
P ( W , λ ′ ( t ) = ∞ ) = exp { − ( W 2 / 1 6 2 − W ) } 4 W π e 32 ( I − 1 / 4 ( 32 ) + I 1 / 4 ( 32 ) ) P(W, \lambda'(t)=\infty)=\frac{\exp \left\{-\left(W^{2} / 16^{2}-W\right)\right\}}{4 \sqrt{W} \pi e^{32}\left(I_{-1 / 4}(32)+I_{1 / 4}(32)\right)} P ( W , λ ′ ( t ) = ∞ ) = 4 W π e 32 ( I − 1/4 ( 32 ) + I 1/4 ( 32 ) ) exp { − ( W 2 /1 6 2 − W ) } 尽管分母中有 W \sqrt W W W ≈ 1 0 2 W\approx 10^2 W ≈ 1 0 2
程序编写
目前写了一组简单的测试程序：
DoubleWell.f90：体系参数和运动方程的求解方法；
Nosé-Poincaré 控温
含时二阶辛算法求解
Expansion.f90：用直接指数平均和累积量展开方法求值；
输出原始和不同阶数展开的估值的期望和方差
基于功的近正态分布讨论误差分析的合理性
数据
每个轨迹在 1 0 − 2 10^{-2} 1 0 − 2 1 0 − 3 10^{-3} 1 0 − 3 1 0 4 10^4 1 0 4
用不同估值方案计算，考察不同估值方案的系统误差；
重复采样 1 0 2 10^2 1 0 2
由于功是正态的，所以二阶展开效果比较好（同时也说明需要更复杂的体系来检验）。
不同方法比较
