如何对 Nosé Hoover 方程进行数值积分?

计算科学 分子动力学
2021-12-14 21:18:42

在 NVT 分子动力学中,Nosé Hoover 恒温器是一种定义扩展系统的方法。我可以完全理解如何推导出 Nosé Hoover 微分方程(Frenkel&Smit 的书)。但是,在实践中,我不知道如何用数值解决它。我尝试了一种类似于velocity Verlet 的方法来整合完整的方程,但它并没有真正起作用。奇怪的是,目前我看到的所有教科书都停留在微分方程上,没有解释如何对其进行数值积分,这才是难点。任何人都可以提出任何可以清楚地解释如何离散化 Nosé Hoover 方程的参考资料,或者任何要查看的书面代码吗?

谢谢!

3个回答

您可以在以下第 4 节中找到 Nosé-Hoover 链的显式伪代码:

Martyna, GJ, Tuckerman, ME, Tobias, DJ 和 Klein, ML (1996)。扩展系统动力学的显式可逆积分器。分子物理学,87(5),1117-1157。取自http://www.tandfonline.com/doi/abs/10.1080/00268979600100761

这是在 MATLAB 中实现的一维谐振子上的 Nosé-Hoover-Langevin 恒温器的简短示例(q是位置,p是动量,ξ是与恒温器相关的自由度,γ是摩擦力,kT是温度乘以玻尔兹曼常数,并且dt是时间步长):

function [q p xi] = nhl(q, p, xi, gamma, kT, dt)

mass = 1.0;

kinetic = @(p) 1/(2*mass) * p^2;

force = @(q) -q;

N = 1;                                  % Number of dimensions

xi = exp(-gamma*dt/2) * xi + sqrt(kT/mass*(1 - exp(-gamma*dt))) * randn;

p = p + dt/2 * force(q);

p = p * exp(-xi*dt/4);
xi = xi + dt/2 * (2*kinetic(p) - N*kT)/mass;
p = p * exp(-xi*dt/4);

q = q + dt/mass * p;

p = p * exp(-xi*dt/4);
xi = xi + dt/2 * (2*kinetic(p) - N*kT)/mass;
p = p * exp(-xi*dt/4);

p = p + dt/2 * force(q);

xi = exp(-gamma*dt/2) * xi + sqrt(kT/mass*(1 - exp(-gamma*dt))) * randn;

您可以找到明确描述的通用 Nosé-Hoover-Langevin 方法:

Leimkuhler, B.、Noorizadeh, E. 和 Theil, F. (2009)。用于分子动力学的温和随机恒温器。统计物理学杂志,135(2),261-277。http://doi.org/10.1007/s10955-009-9734-0

由于这两种方法(NHC 和 NHL)与原始 NH 方法密切相关,因此了解如何在任何这些方法中从 ODE 到实际积分器将引导您了解它是如何为 NH 完成的。

tl;博士

Nose-Hoover 方程通常由它们的哈密顿量定义。哈密​​顿量的两个偏导数定义了一组划分的 ODE。从那里开始,分区 ODE 要么由一阶 ODE 的积分器求解,要么由某种辛方法求解。

使用DifferentialEquations.jl 显示的详细信息

作为一个相当全面的示例,让我们看一下如何处理 DifferentialEquations.jl 的架构。它允许您直接使用哈密顿量来定义动态 ODE 问题。这实际上是在使用哈密顿量的偏导数定义分区 ODE 版本,即由下式给出的位置和动量的运动方程

p˙=Hq
q˙=Hp

来自汉密尔顿的关系。对于方程组,这是梯度,因此您可以从源代码中看到这是如何使用自微分完成的,并在任何语言中类似地定义它。

一旦你有了分区的 ODE 集,你就可以做很多事情。如果将它们组合在一起,则可以将其称为一阶 ODE 的巨型系统。在这种情况下,您可以在该系统上使用任何标准 ODE 求解器(对于 DifferentialEquations.jl ,它是这个列表,或者在 MATLAB 中类似于ode45,或者在 SciPy 中您可以使用LSODA,等等)。

但是您对方程“了解更多结构”,因此理论上可以使用这种结构并更有效地使用方法,对吗?是的,就是这样。在 DifferentialEquations.jl 页面上有一个很大的列表,但我们可以将它们分为两类:Runge-Kutta Nystrom 方法和辛积分器。Runge-Kutta Nystrom 方法是专门为二阶 ODE 开发的 Runge-Kutta 方法的扩展,您在此处拥有的这个结构(位置方程和动量方程)就是这个二阶 ODE 结构。通过删减一些函数评估,在这类特定问题上,这些方法比一阶方法更有效。

辛积分器更有趣。这些积分器满足它们的积分守恒辛二型的性质这通常使人们相信他们节约能源,这是不正确的。然而,对于长期积分,它们将停留在辛流形上,因此不会漂移,从长远来看,可以更好地估计能量和角动量等事物。然而,这些方法有一些与之相关的额外成本(当然,因为它们必须满足更多的属性!)。一件大事是它们是固定的时间步长(有自适应时间步长变体,但最近对该主题的研究并没有表明它们是那么有效,主要是因为自适应时间步长需要继续进行很多额外的工作保持辛)。

回顾

所以让我们回顾一下。您可以通过分析、符号、数字或自微分计算分区 ODE(p˙,q˙): 你选择哪种方法?如果您的问题不太难,只需将其放入标准 ODE 求解器中即可。如果您正在进行长时间积分并且需要“更好”地保留一些额外的属性,例如能量和角动量,请使用一些辛积分器(这里有一个完整的列表,其中一个流行的是 Velocity Verlets,因为它易于实现并减少了力计算,这实际上并不是超级高效,但我会在其他时间进行......)。如果需要它高效且不进行长时间集成,那么您正在寻找一种Runge-Kutta Nystrom 方法。

可用软件

因此,DifferentialEquations.jl 当然如上所示。但是,如果您不使用 Julia,则可以在引文页面上找到每种方法的论文。如果你需要一个辛积分器并且需要自己滚,Velocity Verlet是一个可以先尝试的好东西。除此之外,您可以在 Fortran 中找到其中一些方法,在MATLAB中找到Dormand-Prince Nystrom 方法,在 Python 中找到一些包含 REBOUND的辛积分器,以及在 Mathematica中找到辛积分器,但我不知道还有其他大型一组方法将所有这些都放在一起进行基准测试,因此如果您想要/需要这样做,Julia 是您的选择。

总结性基准和示例

最后,我将介绍一些可以帮助您选择正确方法的基准和示例。这些取自DiffEqTutorials.jlDiffEqBenchmarks.jl,但让我指出一些相关的:

  1. 这是对几个经典问题的辛积分器和标准一阶积分器之间能量误差差异的快速演示
  2. 这个问题显示了在开普勒问题中长时间积分时非辛积分器中发生的“漂移”
  3. 此笔记本在需要高能量精度时使用各种方法对高精度解决方案进行基准测试
  4. 这是另一个高精度基准

虽然这些都不是 Nose-Hoover 方程本身,但这些都是直接从哈密顿量定义的其他问题的基准,并显示了方法类型之间的差异,因此应该可以帮助您从那里推断出您需要什么。

结束快速推荐

如果我不得不在不确切知道您在做什么的情况下提出一项建议,并假设您必须自己实现它(否则……尝试所有这些并进行基准测试!),我会说使用 Velocity Verlet(在分区从哈密顿量定义的 ODE 集)。这里很标准。鉴于我所做的基准测试,symplectic 通常是不必要和过度使用的(人们夸大了 symplectic 的实际含义,在这个类中 Verlets 甚至不是最有效的),但它们很容易实现、标准化,并且可以长期为您提供安全时间整合。

编辑

哎呀,我把它和哈密顿系统的鼻子系统混在一起了。Nose-Hoover 系统是非哈密顿系统,因此标准积分器转换为一阶 ODE 有效,Verlet 方法的一种形式也有效,但标准辛积分器不适用。我们将得到一些实现来“很快”对这个领域进行基准测试。

最简单的方法是使用四阶龙格-库塔积分。