导致每次计算结果不同的原因.docx

1 前言 曾经有很多人问，诸如为什么明明输入文件一样，但是动力学轨迹每次跑出来的都不一样、为什么几何优化结果相差甚远等等，甚至怀疑计算科学是否遵循决定论。 实际上结果的差异是由于运算开始或运算过程中各种形式的数值误差、随机性引起的。这个是很重要却不被多数人重视的问题，它和理论本身、算法、软件环境、硬 件环境等都有密切关联。本文将对这个问题做一些分析探讨。首先先讨论数值算法、软件、硬件因素是如何导致误差（错误）和随机性的。然后再看这些问题对计算化学会产生何等的影响。 2 数值误差（错误）、随机性的根源 2.1 内存因为不少人对内存有误解，所以下面说得多一些。内存在数据读、写的过程中不可避免地会发生一些错误，轻则数据出现异常（比如一个变量为1却成了0），重则死机、重启。和内存错误（包括各种类型错误）相关的主要有以下几个方面： (1)供电质量差。如电压偏离标准值过多、电流不稳。这不仅取决于电源质量，还取决于主板好坏（内存供电电路）。 (2)过热。随着温度升高内存出错的几率总是增加的，但只有高温下增加的幅度才最明显。遥想笔者用i860+RDRAM的日子，不装个小风扇吹着内存的话高负载下20分钟左右就死机。烫手的话应该装上散热片乃至小风扇。 (3)内存质量。一方面是内存设计，即布线、PCB板数、贴片元件等。另一方面取决于内存颗粒制造商的品质，也有RP的关系，比如颗粒来自晶圆的边角还是中心（杂质含量不同）。 (4)制程。一般来说制程越精细可靠性越低，比如22nm制程的比起0.35微米制程的更容易出问题。但不要忽视工艺质量也在逐步改进，新的精细制程的产品出错几率不一定比以往的粗糙制程的要大。 (5)宇宙辐射。这看起来比较微妙，但是当严重的时候，如太阳风暴剧烈时，对电子设备的影响就明显了。这也是为什么在一些宇航领域（接受的辐射更严重） 的芯片仍然用着设计过时、制程古旧的产品，没准儿出现一个小错大家就没命了。机箱内电磁辐射也对内存稳定性有影响。 (6)内存参数。所设的工作频率越高、延迟越低出问题的几率越高，当然也要看内存类型。 (7)内存容量、读写次数：内存容量越大、读写越频繁，相同时间内出错的几率越大。 (8)老化。在越恶劣的环境使用越长时间，出错的可能性越大。要注意以上的因素不是分立、绝对的，而是互相关联的、相对的，不能单独拿出来对比。可见，与内存数据出错相关的因素很多，也很难衡量内存平均多长时间出一次错。一个令我比较认可的、由大规模 统计得到的说法是一般的使用环境下（外界环境正常，计算机运行稳定），平均每年每条内存有3%几率出错（包括各种形式的错误）。实际上这个几率是很小的， 换句话说，等到机子被淘汰了也未必能赶上一次内存出错。所以不要把内存问题过分夸大，比如持续几天的计算，计算结果无法重现由内存错误导致的可能性微乎其 微。但是，对于数据中心，内存出错几率乘上相应内存数目，再考虑数据出错的损失，这就是大问题了。内存的Error Checking and Correcting (ECC)显得十分必要，这是一种数据错误检测和校正的技术，每次在写入数据时同时写入额外的数据位用作校验，在读出数据时也同时把校验位读出，由此根据 一定算法分析读出的数据是否和写入的数据相同，若不同则进行修复。这使得内存数据出错几率大大降低，但也并非能解决所有问题。有些人误认为ECC内存性能 比普通同类型内存要低很多，其实相差甚微。现今服务器内存几乎都带ECC，它比普通内存略贵，多了额外内存颗粒用于储存校验位，而现今服务器主板也都支持 ECC（有些主板则只支持ECC内存）。个人计算没有绝对的必要使用ECC内存。导致内存出错的原因对于CPU/GPU运算出错也基本是一样的，就不累述了。 2.2 处理器类型 对于不同架构的CPU、拥有不同指令集的CPU、乃至GPU，不能期望它们能完全重现同样的结果。尤其是浮点运算的标准可能不同，比如末位的舍入、极小值 的处理。有些指令在不同处理器中执行方式不同，比如先乘一个数再加上一个数，这两条指令在有的处理器中会被融合为结果更精确的“乘加”一条指令。也有一些很特殊的情况，数值错误是由CPU设计bug引起的。老玩家们应该还记得，较早一批的Pentium 60/66Mhz的浮点除法指令有错误。 2.3 网络数据传输 对于集群，尤其是分布式计算，需要考虑网络连接引入的错误。网络传输都会有一定误码率，即在多少位数据中出现一位差错，其大小取决于网络规范、传输 介质、传输距离、外界干扰、网络适配器质量等因素。主流的以太网、infiniband等网络标准本身都带校验帧以保证数据完整、准确性，已很大程度上减 轻了网络传输引发的数据错误。 2.4 随机数生成器 在一些程序中要用到随机数生成器，比如Monte Carlo过程、动力学的初始速度设定。通常随机数是以不同的随机数算法产生的，常用的比如乘同余法。这些算法都需要一个称为“种子”的数，不同的种子将 产生出不同的随机数序列。也就是说，如果程序中利用了随机数生成器，想重现结果就必须用相同的种子。但种子往往不允许由用户设定，也不是固定不变的，而是 根据一些规则产生来保证每次运行结果的随机化，比如利用当前系统的时间，对这种情况用户不要指望结果会重现。 2.5 库文件、操作系统、程序版本 不同的数学库中相同的功能在内部会使用不同的算法或者不同的编写方式，即便是同一个库，在哪怕相近版本中，也可能由于对代码进行的优化、除虫，导致不同 结果。有些库函数的运行结果甚至取决于运行期的情况，比如快速傅里叶变换库FFTW会在启动时检测哪种算法在当前系统下效率最高，然后会在之后一直使用， 这就可能在检测的时候由于正在运行的其它任务对系统CPU资源占用不同，使检测结论不同，影响后续算法。有些程序是通过静态链接产生的，库文件会被合并在程序可执行文件当中，或者有些程序虽然使用动态链接，但是 自带了动态库文件，这样不同系统下结果的重现性还好说。而有些程序需要依赖于当前系统中的库文件，不同系统自带版本不同，差异就容易出现了。若平台不同， 如一个程序的Linux版本与Windows版本，差异会更明显。不要指望不同版本号的程序能获得相同结果，不仅所用的库往往不同，默认参数、算法也经常不同，比如 Gaussian03的几何优化默认用的berny方法，而09版后来默认用GEDIIS，同一个输入文件有时前者能优化成功而后者却不收敛。即便算法和 参数都不变，代码细节也可能改变。 2.6 编译器与编译选项 编译器类型及具体版本号、编译参数对结果影响明显。不同编译器会对代码进行不同形式的优化，有些优化是以牺牲精度为代价的。如果调用了程序语言标准中定义 的标准数学函数，比如开方、求三角函数、求对数等，不同编译器会用不同的内建数学库，显然也会带来差异。不同编译器对于语言标准中没有定义的情况的处理也 不同，比如未初始化的变量，有些编译器一律将之初始化为0，而有些编译器则不这么做，相应内存地址位置目前是什么数值就是什么数值，由于这个数值无法预 测，每次运行时可能造成结果不同。来看一个简单的比较，用Multiwfn 2.02对一个水分子全空间积分电子密度拉普拉斯值的结果。Intel visual fortran与CVF结果有所不同，在不同优化参数下结果也有所不同。IVF120d/O1 -2.603320357104659E-005IVF12O3 -2.603320356662977E-005CVF6.5 debug/release -2.603320356424300E-005 2.7 并行执行时的运算顺序虽然诸如a+b+c=a+(b+c)在数学上是成立的，但是在计算机执行中未必成立。比如在CVF6.5里执行这样的语句a=sin(4D0)b=sin(99D0)c=sin(2D0)write(*,*) a+b+c-(a+(b+c)得 到的结果不是0.000000000000000E+000，而是1.110223024625157E-016，如果结果被乘上一个比较大的数，那么这 个差异就不能被忽略了，可见计算机运算未必满足结合律。在并行运算中，经常要把一部分工作拆成几个线程以分配给多个处理器并行执行，然后将算出来的结果累 加到一起。通常哪个线程先算完，哪个结果就先被累加上。然而实际计算中各个线程完成的先后顺序是不确定的，每次程序执行时都可能不同，因此数据累加的顺序 不一样，这就造成结果出现一定不确定性。另外，有的程序中使用动态均衡负载以获得更高的效率，这导致每个线程执行的内容也是有随机性的，使结果的重现更加 困难。 举个实例，在Gaussian09中用4核计算一个体系单点能，执行三次，最后一次迭代的能量分别为-234.579551993147-234.579551993143-234.579551993144 结果在末位有所差异。如果只用nproc=1来串行计算，则结果每次都能重现。尽管并行时每次的结果差异看似微不足道，以a.u.为单位一般精确到小数 后5、6位足矣，然而这个差异在其它一些任务中会被逐渐放大，甚至可能达到影响定性结论的地步。有些人以为通过增加积分精度、收敛精度之类能减轻结果的随 机性，实际上从原理上讲根本于事无补，根本不是一码事（仅在极个别情况下能有改善）。 2.8 文件格式 不同文件格式的精度，以及格式转换过程中的信息丢失是需要注意的。浮点数据在计算机中通过二进制形式记录，而在磁盘上记录的文件往往是文本格式，这便于 人类阅读和第三方程序处理。在内存中占8字节的双精度变量若想完整以文本方式记录下来起码要用22字节，这不仅太占空间，阅读也不方便，因此往往故意减少 有效数字位数。比如Gaussian的chk文件是二进制文件，浮点变量以二进制形式精确保存，有效数字位数约16。而转换成fch后，浮点数有效数字就 被截为9位，若用unfchk再将之转化回chk，则此时的chk的数据精度比起之前的chk已经损失了。即便是二进制文件也并非都是完整精度形式记录数 据，比如Gromacs的xtc文件也是二进制形式，而变量记录精度可控（xtc_precision参数），默认的记录精度低于单精度浮点变量（这也是 为什么Gromacs轨迹文件比较小的原因），因此不要指望rerun的过程，即重新算轨迹中每一帧的势能的值会与之前动力学过程中输出的一致。有些程序 在导出/导入的过程中还可能会对结构、导数等等进行坐标变换，对单位进行转换等操作，也会引入数值误差。 2.9 操作过程 一些操作过程的细节差异往往不被注意，对最终结果的影响可能却不小。例如对晶体测得的结构加氢，很多程序加氢结果表面看似都一样，位置没有区别，但是查看 文件中记录的坐标，却总有细微差别。再比如模拟一个NMR测得的体系，pdb文件通常会有很多帧，往往每一帧差异很小，似乎取哪帧是随意的，但实际上用于 模拟会得到不同结果。还比如Gaussian中设定不同内存大小，看似和计算结果无关，但实际上Gaussian会自动根据分配的内存容量选择不同的积分 计算、储存方式，也会引入差异。 总结一下，上面的讨论中，数值差异可以归纳为两部分，一种是可重现性差异，也可以叫静态差异，包括硬件配置、编译过程、库文件、文件操作、操作过程、程 序版本，只要运行条件完全一致就可以彻底避免；另一种是不可重现性差异，或者叫动态差异，它取决于运行期的实际条件，有的可以通过一定对策避免，比如把并 行运算改为串行，但有的很难完全消除，尤其是硬件偶发错误。 3 数值误差、随机性对计算化学结果的影响 从前面简单例子看到，误差或者随机性对结果的影响甚微，在默认的输出精度下甚至看不出任何差异。然而对于某些任务情况却不是如此，数值的微小扰动（下面 用扰动这个词代表前述各种因素造成的微小数值差异）可能造成结果有明显定量乃至定性的差异，这是因为在这些任务中扰动的效应会被放大。放大的方式主要可以 归结为两类，一类是数值运算方式，另一类是体系和任务的混沌效应。 3.1 数值运算方式对扰动的放大 一些任务中，由于算法的原因，初始的扰动会在结果中被放大，典型的例子是数值差分。在几何优化、频率计算、超极化率计算等任务中，若理论方法在当前程序中不支持能量解析导数计算，就必须要通过数值差分来近似计算导数。x处一阶导数计算方法为： d1(x)=( E(x+)-E(x-) )/(2) 其中是差分步长。假设最糟的情况，计算E(x+)和E(x-)时都引入了大小为k的扰动，且符号相反而无法被抵消，则此时的一阶导数写为 d1(x)=( E(x+)-E(x-)+2k )/(2)那么由于扰动，造成一阶导数误差为d1(x)-d1(x)=k/。由于必须足够小以保证导数计算准确，比如为0.0001，那么初始扰动k就会被放大10000倍，对结果造成了10000k的影响！ 再来看差分计算二阶导数 d2(x)=( d1(x+)-d1(x-) )/(2) 由于计算d1的时候已经引入了k/的误差，则实际计算出的二阶导数为 d2(x)=( d1(x+)-d1(x-)+2k/ )/(2) 二阶导数的计算误差为d2(x)-d2(x)=k/2，若=0.0001，则初始扰动k就会被放大100000000倍！可见，差分计算n阶导数 时由于初始扰动值k最多可能导致结果产生k/n的误差。可想而知，计算三阶导数和四阶导数时，比如有限场方法纯数值计算一阶和二阶超极化率时，误差已 经被放大到不可忽略的地步了。所以，取得越小，由于误差会越大，结果未必越精确，必须选择最合适的值。通过选择合理的算法、良好的编程习惯，可以使一些任务中误差的放大尽可能小，比如避免大数除小数之类的情况。这就是编程者的事情了。 3.2 体系、任务的混沌性对扰动的放大混沌效应对扰动的放大是特别需要关注的。混沌效应粗俗地说就是捷克一个人从体内释放一股不圣洁的气体可能造 成韩国一场台风，用古话来讲则是“失之毫厘，谬以千里”。混沌效应对数值计算巨大影响的发现，可追溯到1959年洛伦兹用计算机进行气象模拟（其实更早之 前有人在计算天体运行轨道时也多多少少意识到了这点）。他原先模拟时一直使用十几位精度的浮点数保存数据，后来有一次为了更细致地分析某段时间的数据，他 用以前的数据为初值重算了一遍，为了打印的数据美观，这次他将小数精度只保留三位。正是这微小误差的引入，导致模拟出的两个月之后的气象结果和原先大相径 庭！在计算化学中，有两类任务混沌效应很明显，一个是分子动力学，另一个是几何优化。 3.3 扰动对分子动力学的影响分子动力学的混沌效应，尤其是对蛋白质折叠过程的影响被研究得比较多，如PROTEINS,29,417。 动力学过程中初始扰动会随着模拟时间指数型放大，这并不是由于数值计算问题所导致的，而是体系内在性质的体现。若在t=0时对体系坐标引入微扰，则模拟 到t时轨迹与未引入微扰时的差异可写为diff(t)=*exp(t)，其中是依赖于体系的Lyapunov指数。扰动的放大速度相当快，研究发现 对蛋白质结构引入RMSD=1*10-9埃的扰动，经过2ps的模拟就能被放大到RMSD=1埃。假设在模拟中观察到某时刻有一个氨基酸侧链发生了扭 转，如果在此10ps前引入10次RMSD=0.001的不同方式的扰动，并分别跑10次动力学，则只有一次能够重现侧链扭转。可见，动力学过程对初始条 件十分敏感，各种因素带来的微小数值扰动都足矣使动力学轨迹发生很大分歧。更何况实际模拟中扰动并不仅仅是在初始条件中出现，在模拟过程中还可能会不断纳 入。 注意对于蛋白质体系，扰动不是随模拟的进行而无限放大的，因为蛋白的原子受到力场的束缚，且模拟时的温度有限，因此势能面可及的范围是有限的，故初始扰 动造成的差异在迅速放大后会遇到上限。对于柔性体系，由于势能面可及范围会加更大，可以预见差异放大的上限会比刚性体系更高。一些办法可以提高轨迹重现性。提高浮点数精度是重要策略之一，这可以减少扰动量的数量级。有人抱怨 gromacs轨迹重现性不如Amber，这很大程度上是由于gromacs默认是单精度编译，而Amber默认用双精度。不过，提高精度并不解决根本问 题，只不过是让轨迹的歧化慢一些罢了，对于长时间的动力学模拟结果的偏差仍然是无法预测的，除非能让浮点精度无穷高，但这显然是不可能的。若想让轨迹完全 重现，通用的办法是将第一节涉及的所有可能带来扰动的因素都去除。而有些程序也专门考虑了轨迹重现问题，在gromacs里的mdrun有-reprod 选项以帮助轨迹重现。另外还有人弄了个固定精度的分子动力学方法，不仅可以完全重现轨迹，还可以将时间反转从末态得到初始结构，完全体现出分子动力学的决 定性特点。 动力学轨迹的重现究竟有无意义？如果只是出于计算目的，比如一段轨迹文件弄丢了，或者希望以更高的频率输出体系能量，那么重新跑出一次完全一样的轨迹是 有意义的。然而，从理论角度上来讲，轨迹的严格的重现性是毫无意义的。动力学过程就是在体系相空间（坐标+动量空间）中游历的过程，如果跑两次得到了截然 不同的轨迹，只不过是说明在相空间不同位置区域游历罢了，更没有谁对谁错之分。对于一个平衡态的模拟，体系的性质（比如内能、扩散系数）就是动力学过程时间的平均，两次跑出来轨迹相不相 同完全无所谓，只要时间足够长，它们都足以充分遍历相空间，结果会是等价的。对于有限时间的模拟，为了加强某些区域的采样，实际上还有人专门利用动力学轨 迹的混沌效应，给初始速度分配不同的值，分别做动力学模拟并将结果合并。而对于非平衡过程的模拟，比如说蛋白质折叠过程，轨迹（折叠路径）的重现性同样毫无意义。因为蛋白的折叠过 程本来就是混沌的，初始的微小扰动会使折叠路径相差迥异。在现实中每次蛋白的初始条件绝不可能完全一致，因此同样的折叠路径在现实中无法重现，或者说现实 中没有可能发生两次完全相同的折叠过程，那么苛求在多次模拟中重现相同的轨迹有何意义？实际上不仅不应该重现，折叠路径的多样性才是真正感兴趣的，也很适 合通过概率统计来研究，这就需要做很多次动力学模拟，以得到不同的折叠路径并进行分析。正好，由于各种数值因素的介入，哪怕相同的初始条件都会产生不同的 折叠路径，都省得修改输入文件了，直接反复跑多次就行了。注意尽管每次的轨迹十分不同，但最终却一定会折叠到同样的最稳定结构，这并不违背混沌效应，而是 由于每一次折叠过程都会感受到相同的“吸引子”，即趋向成为最稳结构的一种特殊势场。折叠过程可以比喻为站在一个坑坑洼洼的大坑（对应自由能面的形状）旁 边往坑里扔石子，每次石子掉落的轨迹都很不同，但最终都掉到坑底。 这里顺便谈谈一个相关问题，也就是单个动力学模拟的轨迹是否足以说明问题。上面说了，一条折叠轨迹不能反映全部的折叠过程，同样地，研究其它过程，原则上 也需要用多条轨迹来统计地说明问题。比如研究在结合配体后激酶的活性位点打开的过程，在轨迹中发现经过530ps模拟后活性位点打开了，那么就真的能说明 在现实中配体结合后经过约530ps后活性位点就能打开？非也，模拟10次，会得到10个不同的结果，恰好都出现在530ps是不可能的。经过10次模拟 后，假设根据统计结果可说明位点打开的时间主要分布在300ps1.0ns区间内，且最容易在500ps附近打开，则这样的统计性结论比起从单个轨迹中 下的结论有意义、可信得多。混沌效应导致动力学过程中各种“事件”发生的时刻的波动范围极大，或者说每种事件在不同的时间范围都有一定出现概率，只分析单 一轨迹就造成了一叶障目，结论的不可靠性会较大。然而目前绝大多数模拟都是基于单一轨迹的，很少有人对同一问题做多轨迹统计分析，主要是多轨迹的计算量太 大或者研究的内容不需要那么深入详细，其次是很多人没意识到这个问题，还有少数人是因为投机取巧。比如论证一个抑制剂的作用，可能模拟反复跑了7次都无法 论证抑制剂管用，当跑到第8次时庆幸地发现能说明抑制剂管用，于是他只用这一次的轨迹说事，另外跑的7次虽然没成功却都避而不谈，可想而知这抑制剂的实际 功效如何了. 3.4 扰动对几何优化的影响相较分子动力学的混沌效应，几何优化的混沌效应被关注的相对较少，有兴趣者可参看J Comput Aided Mol Des,22,39，讨论数值误差对分子力场优化结果的影响。之所以几何优化也有明显混沌效应，不妨将几何优化看作是0K下的分子动力学，尤其对于最陡下降法优化，这个 比喻最为恰当。最陡下降法优化过程类似把一个玻璃珠轻轻放到陡峭的山上某斜坡处，让它自动滚落。虽然0K下没有动能（准确来说尚有振动能），通常不会导致 引入扰动后几何优化过程的轨迹像动力学轨迹那样随步数增长分歧得那么快，但是，由于没有了动能以越过势垒，体系将没机会感受到全局的吸引子（全局最稳结 构），而只能感受到离初始位置最近的局部吸引子（势能面局部极小点），这导致了几何优化的最终结果的重现性可能比动力学更差。对于小体系还好，势能面结构 比较简单、光滑，然而极小点、鞍点的数目是随着原子数目增加呈指数型增长的，对于大体系，势能面的复杂程度是难以想象的，微小的扰动可能就会使结构收敛到 很不同的极小点。之所以平时在量子化学的几何优化中优化的结果重现性都很好，没有显示出太大不可重现性，主要是因为量化处理 的体系一般都比较小，通常少于50个原子，所需优化步数一般也较少，通常少于100步。不仅误差来不及充分放大，而且由于相对简单的势能面的形状以及它起 到的束缚作用，误差也不容易放大。然而因数值误差带来的分歧偶尔仍可能显著，记得有人在一台机子上某个体系优化几十步成功了，换了台机子却不收敛了。重现性好坏与初猜坐标关系很大，在山脚下释放一个珠子，和在山尖上释放一个珠子，明显珠子最终位置在后者的 情况中受初始位置影响更大。所以优化时如果能给出一个好的初猜，比如在极小点附近的二次区域内（势能面可以近似用二次函数描述的区域），那么优化几乎一定 收敛到那个极小点，数值误差因素不会导致收敛到别处去。然而如果在很接近过渡态的位置作为初猜，一开始的坐标或者计算受力稍微有点扰动，很可能结构就会收