董亚琴　王圣军

摘要：本文以简谐振子为例讨论在计算机模拟中保守系统的能量变化。模拟和解析分析表明常微分方程通用解法使能量随着时间指数地变化，而Verlet算法只带来振子能量的涨落，能量不会单调变化。这些结果为理解分子动力学模拟算法的特性提供了直观的参考。

关键词：分子动力学；误差；能量

中图分类号：O411 文献标志码：A 文章编号：1674-9324（2018）50-0193-02

在分子动力学教学中一个重要的内容是选择求微分方程的数值解法。尽管存在多种通用解法，能够求解各种形式的常微分方程，但是它们在求解力学系统的运动方程时不能保证能量守恒。分子動力学模拟需要避免这个问题。虽然能量误差在分子动力学中非常重要，但是各种教材只有简单的描述。例如参考文献[1]指出：欧拉方法的时间步长太大时，简谐振子的相空间轨道不是椭圆，能量被粒子吸收。参考文献[2]指出：用二阶龙格库塔方法模拟简谐振子的运动时其相空间轨道收缩。参考文献[3]提到：模拟简谐振子运动时步长取得太大得到的结果偏离也大，甚至连能量都不守恒。这些说法都不容易理解，甚至不够准确。为帮助学生深入理解相关知识，本文讨论简谐振子模拟中通用解法的能量误差。

模拟得到的能量时间关系如图2。它不引起能量的单调变化，只带来能量涨落，步长越小涨落的范围越小。

本文的结果表明欧拉方法、泰勒方法、二阶龙格库塔方法导致谐振子能量随着时间单调递增，相空间轨道螺旋增大。用四阶龙格库塔方法模拟时能量随时间单调递减，减小步长或者改变参数不能改变能量的变化趋势。使用Verlet算法能量只发生涨落。这些结果详细说明了通用算法在处理牛顿运动方程时的误差特性，为解释分子动力学中算法的选择提供了详细的实例和分析，是教材的有益补充。

参考文献：

[1]D.W.Heermann.理论物理学中的计算机模拟方法[M].北京：北京大学出版社，1996：24-25，5-6.

[2]顾昌鑫.计算物理学[M].上海：复旦大学出版社，2010：396-397.

[3]马文淦.计算物理学[M].北京：科学出版社，2005：92-93.

[4]S.E.Koonin.计算物理学[M].北京：高等教育出版社，1993：19-20.

[5]N.J.Giordano，H.Nakanishi.计算物理[M].北京：清华大学出版社，2007：273.