N年卜洗澡
【物理解题一般程序】
【物理解题方法】
整体法、假设法、极限法、逆向思维法、物理模型法、等效法、物理图像法,等. 【典型习题】
1.质量为 2kg 的物体静止在足够大的水平地面上,物体与地面间的动摩
擦因数为 0.2,最大静摩擦力与滑动摩擦力大小视为相等。从 t = 0时 刻开始, 物体受到方向不变、 大小呈周期性变化的水平拉力 F 的作用, F 随时间 t 的变化规律如图所示。重力加速度 g 取 10m/s2, 则物体在 t = 0至 t = 12s这段时间的位移大小为 ( )
A . 18m B . 54m C . 72m D . 198m 2.一个静止的质点,在 0 ~ 5s 时间内受到力 F 的作用,力的方向始终 在同一直线上, 力 F 随时间 t 的变化图线如图所示。 则质点在 ( ) A .第 2s 末速度方向改变 B .第 2 s末加速度为零 C .第 4 s末运动速度为零 D .第 4 s末回到原出发点
3. 将一个物体以某一速度从地面竖直向上抛出, 设物体在运动过程中所受空气阻力大小不变, 则物体 A .刚抛出时的速度最大 B .在最高点的加速度为零
C .上升时间大于下落时间 D .上升时的加速度等于下落时的加速度 4.如图,质量为 M 的楔形物 A 静置在水平地面上, 其斜面的倾角为 θ.斜面
上有一质量为 m 的小物块 B , B 与斜面之间存在摩擦. 用恒力 F 沿斜面向 上拉 B , 使之匀速上滑. 在 B 运动的过程中, 楔形物块 A 始终保持静止. 关 于相互间作用力的描述正确的有 ( ) A . B 给 A 的作用力大小为 mg – F B . B 给 A 摩擦力大小为 F
C .地面受到的摩擦力大小为 F cos θ D .地面受到的压力大小为 Mg + mg – F sin θ 5.如图所示, ad 、 bd 、 cd 是竖直面内的三根固定的光滑细杆, a 、 b 、 c 、 d 位于同一圆周上, a 点为圆周上最高点, d 点为圆周上最低点。每根杆上都 套有一个小圆环,三个圆环分别从 a 、 b 、 c 处由静止释放,用 t 1 、 t 2 、 t 3依次表示各环到达 d 点所用的时间,则 ( ) A . t 1
6.如图所示,光滑水平面上有质量分别为 m 1和 m 2的甲、乙两木块,两木块中
间用一原长为 L 、 劲度系数为 k 的轻质弹簧连接起来, 现用一水平力 F 向左
推木块乙,当两木块一起匀加速运动时,两木块之间的距离是 ( ) A . L + m 2F /[(m 1 + m 2) k ] B . L – m 1F /[(m 1 + m 2) k ] C . L – m 1F /(m 2k ) D . L + m 2F /(m 2k )
7.如图所示,物体 ABC 放在光滑水平面上用细线 ab 连接,力 F 作
用在 A 上,使三物体在水平面上运动,若在 B 上放一小物体 D , D 随 B 一起运动,且原来的拉力 F 保持不变,那么加上物体 D
后两绳中拉力的变化是 ( ) A . T a 增大 B . T b 增大 C . T a 变小 D . T b 不变
8.在水平的足够长的固定木板上,一小物块以某一初速度开始滑动,经一段时间 t 后停止.现将该木 板改置成倾角为 45°的斜面,让小物块以相同的初速度沿木板上滑.若小物块与木板之间的动摩擦 因数为 μ.则小物块上滑到最高位置所需时间与 t 之比为 ( ) A .
μ
μ+12 B .
μ
μ
21+ C .
μ
μ
+2 D .
μ
μ21+
9.某人在地面上用弹簧秤称得体重为 490N 。他将弹簧秤移至电梯内称其体重, t 0至 t 3时间段内,弹簧秤的示数如图所示,电梯运行的 v – t图可能是(取电 梯向上运动的方向为正) ( )
10.航模兴趣小组设计出一架遥控飞行器,其质量 m =2㎏,动力系统提供的恒定升力 F =28 N。试飞
时,飞行器从地面由静止开始竖直上升。设飞行器飞行时所受的阻力大小不变, g 取 10m/s2。 ⑴ 第一次试飞,飞行器飞行 t 1 = 8 s 时到达高度 H = 64 m。求飞行器所阻力 f 的大小;
⑵ 第二次试飞,飞行器飞行 t 2 = 6 s 时遥控器出现故障,飞行器立即失去升力。求飞行器能达到 的最大高度 h ;
⑶ 为了使飞行器不致坠落到地面,求飞行器从开始下落到恢复升力的最长时间 t 3。
a
c
d
甲
F
乙
11.如图所示,一质量 M = 0.2kg的长木板静止在光滑的水平地面上,另一质量 m = 0.2kg的小滑块, 以 v 0 = 1.2m/s的速度从长木板的左端滑上长木板。已知小滑块与长木板间的动摩擦因数 μ1 = 0.4, g = 10m/s2,求:
⑴ 经过多少时间小滑块与长木板速度相等?
⑵ 从小滑块滑上长木板, 到小滑块与长木板相对静止,
小滑块运动的距离为多少?(滑块始终没有滑离长木板)
12.放在水平地面上的一物块,受到方向不变的水平推力 F 的作用,力 F 的大小与时间 t 的关系和物 块速度 v 与时间 t 的关系如图所示,重力加速度 g = 10m/s2,求:
⑴ 物块在运动过程中受到的滑动摩擦力大小;
⑵ 物块在 3— 6s 中的加速度大小;
⑶ 物块质量是多少?
13.如图所示,质量为 m = lkg的小球穿在斜杆上,斜杆与水平方向的夹角为 θ= 30°,球恰好能在杆 上匀速滑动.若球受到一大小为 F = 20N 的水平推力作用,可使小球沿杆向上加速滑动(g 取 10m/s2) ,求:
⑴ 小球与斜杆间的动摩擦因数 μ的大小;
⑵ 小球沿杆向上加速滑动的加速度大小. 14.卡车拖挂着相同质量的车厢,在水平直道上以 v 0 = 12 m/s的速度匀速行驶,其所受阻力可视为与 车重成正比,与速度无关。某时刻,车厢脱落,并以大小为 a = 2 m/s2的加速度减速滑行。在车 厢脱落 t = 3s后,司机才发觉并紧急刹车,刹车时阻力为正常行驶时的 3倍。假设刹车前牵引力 不变,求卡车和车厢都停下后两者之间的距离。
15.一小圆盘静止在桌布上,位于一方桌的水平桌面的中央。桌布的一边与桌的 AB 边重合,如图。 已知盘与桌布间的动摩擦因数为 μ1,盘与桌面间的动摩擦因数为 μ2。现突然以恒定加速度 a 将桌 布抽离桌面,加速度方向是水平的且垂直于 AB 边。若圆盘最后未从桌面掉下,则加速度 a 满足 的条件是什么?(以 g 表示重力加速度)
16.一水平的浅色长传送带上放置一煤块(可视为质点) ,煤块与传送之间的动摩擦因数为 μ。初始时 , 传 送带与煤块都是静止的,现让传送带以恒定的加速度 a 0开始运动,当其速度到达 v 0后,便以此速度 做匀速运动。经过一段时间 , 煤块在传送带上留下了一段黑色痕迹后,相对于传送带不再滑动,求此 黑色痕迹的长度。
答案:
1. B 2. C 3. A 4. CD 5. D 6. B 7. A 8. A 9. A 10.⑴ 第一次飞行中,设加速度为 a 1,匀加速运动 H = a 1t 12/2
由牛顿第二定律 F – mg – f = ma 1 解得 f = 4N
⑵ 第二次飞行中,设失去升力时的速度为 v 1,上升的高度为 s 1
匀加速运动 s 1 = a 1t 22/2,设失去升力后的速度为 a 2,上升的高度为 s 2 由牛顿第二定律 mg + f = ma 2,又 v 1 = a 1t 2、 s 2 = v 12/2a 2, 解得 h = s 1 + s 2 = 42m
⑶ 设失去升力下降阶段加速度为 a 3;恢复升力后加速度为 a 4,恢复升力时速度 为 v 3, 由牛顿第二定律 mg – f = ma 3、 F + f – mg = ma 4, 且 v 32/2a 3 + v32/2a 4 = h ,
v 3 = a 3t 3,解得 t 3
=
2
s (或 2.1s ) 11.⑴ 分析 m 、 M 的受力,由牛顿第二定律有 a m = μ1g = 4m/s2、 a M = μ1g = 4m/s2、
设经过时间 t 两者速度相同,则 v m = v 0 – a m t 、 v M = a M t 当 v m = vM 时,代入数据,联解可得 t = 0.15s
⑵ 小滑块做匀减速运动.初速度 v m = v 0 – a m t = 0.6m/s, S = (v t 2 – v m 2)/2a = 0.135m。 12.⑴ 在 6— 9s 内物体匀速运动,则物体所受滑动摩擦力 f = F = 4N;
⑵ 由于 v – t 图象的斜率就是物体的加速度 a = 2m/s2; ⑶ 在 3— 6s 内 F – f = ma 得 m = 1kg。
13.⑴ 对小球受力分析,由平衡条件可知,平行于杆方向:mg sin θ = f 1, y 轴方向:N 1=mg cos θ
f 1 = μN1 解得小球与斜杆间的动摩擦因数 μ = tan30°=
3
3。 ⑵ 水平推力作用后,由牛顿第二定律:F cos θ– mg sin θ– f 2 = ma 、 f 2 = μN2 = μ(F sin θ+mg cos θ)
解得小球沿杆向上加速滑动的加速度:a =
3
3
20– 10=1.55m/s2. 14.设卡车的质量为 M ,车所受阻力与车重之比为 μ;刹车前卡车牵引力的大小为 F ,卡车刹车前后
加速度的大小分别为 a 1和 a 2。重力加速度大小为 g 。由牛顿第二定律有 F – 2μMg = 0、 F – μMg = Ma 1、 μMg = Ma 、 3μMg = Ma 2
设车厢脱落后, t = 3s内卡车行驶的路程为 s 1,末速度为 v 1,根据运动学公式有 s 1 = v 0t + a 1t 2/2、 v 1 = v 0 + a 1t 、 v 12 = 2a 2s 2 式中 s 2是卡车在刹车后减速行驶的路程 设车厢脱落后滑行的路程为 s ,有 v 02 = 2as 卡车和车厢都停下来后相距 Δs = s 1 + s 2 – s
联立解得 Δs = 34320+-
a v v 0t + 3
2
at 2,带入题给数据得 Δs = 36m。 15.设圆盘的质量为 m ,桌长为 l ,在桌布从圆盘上抽出的过程中,盘的加速度为 a 1,有 a 1 = μ1g
桌布抽出后,盘在桌面上作匀减速运动,以 a 2表示加速度的大小,有 a 2 = μ2g
设盘刚离开桌布时的速度为 v 1,移动的距离为 x 1,离开桌布后在桌面上再运动距离 x 2后便停下,
有 v 12 = 2a 1x 1、 v 22 = 2a 2x 2 盘没有从桌面上掉下的条件是 x 2 ≤ l /2 – x 1
设桌布从盘下抽出所经历时间为 t ,在这段时间内桌布移动的距离为 x ,有 x = at 2/2、 x 1 = a1t 2/2
而 x = l /2 + x 1 由以上各式解得 a ≥ (μ1 + 2μ2) μ1g /μ2。
16. 根据“传送带上有黑色痕迹”可知, 煤块与传送带之间发生了相对滑动, 煤块的加速度 a 小于传送带 的加速度 a 0,根据牛顿第二定律可得 a = μg。设经历时间 t ,传送带由静止开始加速到速度等于 v 0,煤 块则由静止加速到 v ,有 v 0 = a0t 、 v = at ;由于 a < a="" 0,故="" v="">< v0,煤块继续受到滑动摩擦力的作用。再="" 经过时间="" t="" '="" ,煤块的速度由="" v="" 增加到="" v="" 0,有="" v="" 0="v" +="" at="" '="" ,此后,煤块与传送带运动速度相同,相对于="" 传送带不再滑动,不再产生新的痕迹。设在煤块的速度从="" 0增加到="" v="" 0的整个过程中,传送带和煤块移="" 动的距离分别为="" s="" 0和="" s="" ,有="" s="" 0="a" 0t="" 2/2="" +="" v="" 0t="" ′="" ,="" s="v" 02/2a="" ,传送带上留下的黑色痕迹的长度="" l="s" 0="" –="" s="" ,由="" 以上各式得="" l="v" 02(a="" 0="" –="" μg="" )/(2μa="" 0g="" )="">
分子动力学方法
分子动力学方法 一、引言
计算机模拟中的另一类确定性模拟方法，即统计物理中的所谓合于动力学方法 （Molecular
Dynamics Method）。这种方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要
建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各
个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性，
从而得到系统的宏观性质。在这样的处理过程中我们可以看出：MD方法中不存在任何随机因素。在MD方法处理过程中方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系
中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者
拉格朗日量来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现Boltzman的统计力学途径。
这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，讨
算量大，占内存也多、本节将介绍分子动力学方法及其应用。
原则上，MD方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，
也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分
子，也可以是其他的微观粒子。
实际上，MD模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间的限制；
另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数日趋于无穷时）的
性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减
小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响
体系的某些性质。
对于MD方法，向然的系综是微正则系综，这时能量是运动常量。然而，当我们想要研究温度和
（或）压力是运动常量的系统时，系统不再是封闭的。例如当温度为常量的系统可以认为系统是放置在
一个热俗中。当然，在MD方法中我们只是在想像中将系统放入热浴中。实际上，在模拟计算中具体
所采取的做法是对一些自由度加以约束。例如在恒温体系的情况下，体系的平均动能是一个不变量。这
时我们可以设计一个算法，使平均动能被约束在一个给定值上。由于这个约束，我们并不是在真正处理
一个正则系综，而实际上仅仅是复制了这个系综的位形部分。只要这一约束不破坏从一个状态到另一个
状态的马尔科夫特性，这种做法就是正确的。不过其动力学性质可能会受到这一约束的影响。
自五十年代中期开始，MD方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的
重要方法。应用MD方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，
以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。
二、分子运动方程的数值求解
采用MD方法时，必须对一组分于运动微分方程做数值求解。从计算数学的角度来看，这个求解
是一个初值问题。实际上计算数学为了求解这种问题己经发展了许多的算法，但并不是所有的这些算法
都可以用来解决物理问题。下面我们先以一个一维谐振子为例，来看一下如何用计算机数值计算方法求
解初值问题。一维谐振子的经典哈密顿量为
（2.1）
这里的哈密顿量（即能量）为守恒量。假定初始条件为x(p)、p(0),则它的哈密顿方程是对时间的一阶微分方程
（2.2）
现在我们要用数值积分方法计算在相空间中的运动轨迹(X(t)、p(t)) 。我们采用有限差分法，将微分方程变为有限差分方程，以便在计算机上做数值求解，并得到空间坐标和动量随时间的演化关系。首先，
1
我们取差分计算的时间步长为h，采用我们有限差分法中的一阶微分形式的向前差商表示，即直接运用
展开到h的一阶泰勒展开公式
即
（2.3）
则微分方程（2.2）可以被改写为差分形
（2.4）
（2.5）
将上面两个公式整理后，我们得到解微分方程(2.2)的欧拉(Euler)算法
（2.6）
（2.7）
这是x(t)、p(t)的一组递推公式、有了初始条件x(0)、p(0)，就可以一步一步地使用前一时刻的坐标、动量值确定下一时刻的坐标、动量值。这个方法是一步法的典型例子。
由于在实际数值计算时h的大小是有限的，因而在上述算法中微分被离散化为差分形式来计算时总
是有误差的。可以证明一步法的局部离散化误差与总体误差是相等的，都为O(h
2)的量级。在实际应用
中，适当地选择h的大小是十分重要的。h取得太大，得到的结果偏离也大，甚至于连能量都不守恒：
h取得太小，有可能结果仍然不够好。这就要求我们改进计算方法，进一步考虑二步法。
实际上泰勒展开式的一般形式
(2.8)
n+1其中O(h)表示误差的数量级。前面叙述的欧拉算法就是取n=1。现在考虑公式（2.8）中直到含h的
二次项的展开（即取n＝2），则得到
(2.9)
(2.10)
将上面两式相加、减得到含二阶和一阶导数的公式
(2.11)
(2.12)
2令f(t)＝x(t)，利用牛顿第二定律公式dxF(t),m，公式(2.11)写为坐标的递推公式 2dt
2
(2.13)
公式(2.12)写为计算动量的公式得到
(2.14)
这样我们就推导出了一个比(2.6)和(2.7)更精确的递推公式。这是二步法的一种，称为Verle方法。还有
其他一些二步法，如龙格一库塔（Runge－Kutta）方法等。
当然我们还可以建立更高阶的多步算法，然而大部分更高阶的方法所需要的内存比一步法和二步法
所需要的大得多，并且有些更高阶的方法还需要用迭代来解出隐式给定的变量，内存的需求量就更大。
并且当今的计算机都仅仅只有有限的内存，因而并不是所有的高阶算法都适用于物理系统的计算机计
算。
在实际数值计算中，我们必须特别注意舍入误差和稳定性问题。为了减少舍入误差，我们可以采用
高精度计算,并且要避免相近大小的数相消，以及数量级相差很大的两个数相加和注意运算顺序。
三、分子动力学模拟的基本步骤
在计算机上对分子系统的MD模拟的实际步骤可以划分为四步；首先是设定模拟所采用的模型：
第二，给定初始条件；第三，趋于平衡的计算过程；最后是宏观物理量的计算。下面就这四个步骤分别
做简单介绍。
1、模拟模型的设定
设定模型是分子动力学模拟的第一步工作。例如在一个分子系统中，假定两个分子间的相互作用势
为硬球势，其势函数表示为
实际上，更常用的是图（3.1）Lennard-Jones型势。它的势函数表示为
（3.1）
其中-ε是位势的最小值（ ε可以确定能量的单位），这个最小值出现在距离 r等于21/6ζ的地方（ζ可以
确定为长度的单位）
图3.1 Lennard-Jones势
模型确定后，根据经典物理学的规律我们就可以知道在系综模拟中的守恒量。例如对在微正则系综的模
拟中能量、动量和角动量均为守恒量。在此系综中它们分别表示为
（3.2）
3
（3.3）
（3.4）
其中p=mr。由于我们只限于研究大块物质在给定密度卜的性质，所以必须引进一个叫做分子动力学元ii
胞的体积元，以维持一个恒定的密度。对气体和液体，如果所占体积足够大，并且系统处于热平衡状态
的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，元胞的形状是有影响的。为了计算简
便，对于气体和液体，我们取一个立方形的体积为MD元胞。设MD元胞的线度大小为L，则其体积
3为L。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱的表面的
粒子应当被反射回到元胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的任何一种
性质都会有重大的影响。为了减小引入的表面效应，我们采用周期性边界条件。采用这种边界条件，我
们就可以消除引入的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。实际上，这里我
们做了一个假定，即让这个小体积元胞镶嵌在一个无穷大的大块物质之中。周期性边界条件的数学表示
形式为
（3.5） 其中A为任意的可观测量，n、n、n。为任意整数。这个边界条件就是命令基本MD元胞完全等同地123
重复无穷多次。具体在实现该边界条件时是这样操作的：当有一个粒子穿过基本MD元胞的六方体表面时，就让这个粒子以相同的速度穿过此表面对面的表面重新进入该MD元胞内。
在分子动力学模拟中考虑粒子间的相互作用时，通常采用最小像力约定。这个约定是在由无穷重复
的MD基本元胞中，一个粒子只同它所在的基本元胞内的另外N-1个（设在此元胞内有N个粒子）中的每个粒子或其最邻近影像粒子发生相互作用。如果r处的粒子i同r处的粒子j之间的距离为 ij
（对一切的n） (3.6)
实际上这个约定就是通过满足不等式条件rc
很大，使得距离大于L/2的粒子的相互作用可以忽略，以避免有限尺寸效应。采用最小像力约定使得在
截断处粒子的受力有一个δ-函数的奇异性，这会给模拟计算带来误差。
2、给定初始条件
MD模拟进人对系统微分方程组做数值求解的过程时，需要知道粒子的初始位置和速度的数值。不
同的算法要求不同的初始条件。例如，Verlet方法需要两组坐标来启动计算；一组是零时刻的坐标，另
一组是前进一个时间步长时的坐标，或者是一组零时刻的速度值。但是，一般来说系统的初始条件都是
不可能知道的。表面上看这是一个难题。实际上，精确选择待求系统的初始条件是没有什么意义的，因
为模拟时间足够长时，系统就会忘掉初始条件。但是初始条件的合理选择将可以加快系统趋于平衡。常
用的初始条件可以选择为：（1）令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机
抽样得到。（2）令初始位置随机地偏离差分划分网格的格子，初始速度为零。（3）令初始位置随机地偏离差分划分网格的格子，初始速度从玻尔兹曼分布随机抽样得到。
3．趋于平衡
按照上面给出的运动方程、边界条件和初始条件，就可以进行分子动力学模拟计算。但是，这样计
算出的系统不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为了使系统达到平衡，
模拟中需要一个趋衡过程。在这个过程中，增加或从系统中移出能量，直到系统具有所要求的能量。然
后，再对运动方程中的时间向前积分若于步，使系统持续给出确定能量值。我们称；这时系统已经达到
平衡态。这段达到平衡所需的时间称为弛豫时间。在MD模拟中，时间步长h的大小选择是十分重要的。它决定了模拟所需要的时间。为了减小误差，步长h必须取得小一些；但是取得太小，系统模拟的
弛豫时间就很长。这里需要积累一定的模拟经验，选择适当的时间步长h。例如，对一个具有几百个氩( Ar)分于的体系，如果采用Lennard-Jones位势，我们发现取h为10-2量级，就可以得到好的相图. 这里
-14选择的h是没有量纲的，实际上这样选择的h对应的时间在10秒的量级。如果模拟1000步，系统达
-11到平衡态，驰豫时间只有10秒。
4
4．宏观物理量的计算
实际计算宏观物理量往往是在MD模拟的最后阶段进行的。它是沿着相空间轨迹求平均来计算得
(N)到的。例如对于一个宏观物理量A，它的测量值应当为平均值。如果已知初始位置和动量为{r(0)}A
(N)和{p(0)}（上标N表示系综N个粒子的对应坐标和动量参数），选择某种MD算法求解具有初值问题
(N)(N)的运动方程，便得到相空间轨迹{r(t)}，{p(t)}对轨迹平均的宏观物理量A的表示为
(3.7)
如果宏观物理量为动能，它的平均为
(3.8)
由于在模拟过程中计算出的动能值是在不连续的路径上的值，因此公式（3.8）可以表示为在时间的各个间断点μ上计算动能的平均值
（3.9）
在MD模拟过程中，温度是需要加以监测的物理量，特别是在模拟的起始阶段。根据能量均分定理，
我们可以从平均动能值计算得到温度值：
（3.10）
其中d为每个粒子的自由度，如果不考虑系统所受的约束，则d= 3。系统内部的位形能量的轨道平均值为
（3.11）
假定位势在r
处被截断，那么上式计算出的势能以及由此得到的总能量就包含有误差。为了对此偏差C
作出修正，我们采用对关联函数来表示位能
（3.12）
式中的g(r)就是对关联函数，它是描述与时间无关的粒子间关联性的量度。g(r) 的物理意义是当原点r=0处有一个粒子时，在空间位置r的点周围的体积元中单位体积内发现另一个粒子的几率。若n(r)为距离原点r到r+Δr之间的平均粒子数，则
（3.13）
在MD模拟过程中，所有的距离已经在力的计算中得到，因而很容易计算对关联函数的值。图3.2为由计算机模拟得到的两组不同参数下的对关联函数的例子。由于位势的截断，对关联函数仅对r＜L/2以C
下的距离有意义。在公式（3.11）中，所有的位能都加到截断距离为止，尾部修正可以取为
（3.14）
压强可以通过计算在面积元dA的法线方向上净动量转移的时间平均值来得到，也可以利用含对关联函
数的维里状态方程计算。该维里状态方程可以写为
5
（3.15）
至于势能的计算，我们可以把积分划分为两项，一项是由相互作用力程之内的贡献引起的，一项是对位
势截断的改正项：
（3.16）
其中长程改正项为：
（3.17）
下面我们将讨论具体如何进行MD模拟：
图6.2 由计算机模拟得到的两组不同参数下的对关联函数 4、平衡态分子动力学模拟
在经典MD模拟方法的应用当中，存在着对两种系统状态的MD模拟。一种是对平衡态的MD模
拟，另一类是对非平衡态的MD模拟。对平衡态系综MD模拟又可以分为如下类型：微正则系综的MD
（NVE）模拟，正则系综的MD（NVT）模拟，等温等压系综MD（NPT）模拟和等焓等压系综MD（NPH）
模拟等。下面我们仅对平衡态的MD方法中前两类模拟做简单的介绍。
1）、微正则系综的MD模拟
在进行对微正则系综的MD模拟时，首先我们要确定所采用的相互作用模型。我们假定一个孤立
的多粒子体系，其粒子间的相互作用位势是球对称的，则其哈密顿量可以写为
（4.1）
其中r
是第i个粒子与第j个粒子之间的距离。在这个微正则系综中，由于这个系统的哈密顿量中不显ij
式地出现时间关联，因而系统的能量是个守恒量。系统的体积和粒子数也是不变的。此外；由于整个系
统并未运动，所以整个系统的总动量P恒等于零。这就是系统受到的四个约束。
由该系统的哈密顿量可以推导出牛顿方程形式的运动方程组
(4.2)
要用数值求解的方法解出（4.2）微分方程组，类似于本章第二节中介绍的Veriet方法方程组（4.2）的求解变成求解方程组：
(4.3)
6
该方程组反映出：从前面t和t-h时刻这两步的空间坐标位置及t时刻的作用力，就可以算出下一步 t+h时刻的坐标位置。下面为了将（4.3）式写成更简洁的形式，令
(4.4)
则从（4.3）式可以得到如下差分方程组的形式
(4.5)
(0)(1)(2)(3)如果已知一组初始空间位置{{r},{r}{r}{r}，则通过求解方程组（4.5）一步步得到,,……。iiii由空间坐标又可以算出粒子的运动速度为
(4.6)
这里在第 n＋ l步算出的速度是前一时刻；即第 n步的速度、因而动能的计算比势能的计算落后一步。
根据上述原理我们可以将粒子数恒定、体积恒定、能量恒定的微正则系综（NVE）的MD模拟步骤设计如下：
（1）给定初始空间位置(0)(1){{r},{r}，(i=1,2,3,……N) ii
(n)(n) （2）计算在第 n步时粒子所受的力F,F(t){F}：。 iiin
(n，1)(n)(n,1)(n)2(n，1) （3）利用公式r,2r,r，Fh/m{r}计算在第n+l步时所有粒子所处的空间位置 iiiii
(n)(n，1)(n,1) （4）计算第n步的速度：v,(r,r)/2h iii
（5）返回到步骤（2），开始下一步的模拟计算。
如前所述，用上述形式的Verlet算法，动能的计算比势能的计算落后一步。此外，这种算法不是自
(0)启动的。要真正求出微分方程组（4.2）的解，除了需要给出初始空间位置{r}外，还要求另外给出一i
(1)组空间位置{r}。实际，有时候采用改进后的计算方法可能更方便；即把N个粒于的初始位置放在网i
(1)格的格点上，然后加以扰动。如果初始条件是空间位置和速度，则采用下面的公式来计算空间位置{r} i
（4.7）
然后再按上述模拟步骤进行计算。
verlet算法的速度变型形式将会使其数值计算的稳定性得到加强。下面我们就此做简单介绍。令
（4.8）
则公式（4.5）写为
（4.9）
上式在数学上与（4.5）式是等价的，并称为相加形式。由此Verlet算法的速度形式的模拟步骤可以表
述为
（l）给定初始空间位置(1){r}，(i=1,2,3,……N) i
7
(1) (2） 给定初始速度 {v}i
(n，1)(n)(n)(n)2 （3）利用公式：r,r，hv，Fh/2m，计算在第n＋1步时所有粒子所处的空间位置iiii
(n，1){r}。 i
(n，1)(n，1)(n)(n，1)(n) （4）计算在第n+1步时所有粒子的速度{v}v,v，h(F，F)/2m： iiiii
（5）返回到步骤（3），开始第n+2步的模拟计算。
Verlet速度形式的算法比前一种算法好些。它不仅可以在计算中得到同一时间步上的空间位置和速
度，井且数值计算的稳定性也提高了。
一般情况下，对于给定能量的系统不可能给出精确的初始条件。这时需要先给出一个合理的初始条
件，然后在模拟过程中逐渐调节系统能量达到给定值。其步骤为：首先将运动方程组解出若干步的结果；
然后计算出动能和位能；假如总能量不等于给定恒定值，则通过对速度的调整来实现能量守恒。也就是
将速度乘以一个标度（scaling）回子，该因子一般取为
（4.10）
然后再回到第一步，对下一时刻的运动方程求解。反复进行上面的过程，直到系统达到平衡。 这样的模拟过程也称为平衡化阶段。
采用对速度标度的办法，可以使速度发生很大变化。为了消除可能带来的效应，必须要有足够的时
间让系统再次建立平衡。在到达趋衡阶段以后，必须检验粒子的速度分布是否符合麦克斯韦- 波尔兹曼分布。
2）、 正则系综的MD模拟
在统计物理中的正则系综模拟是针对一个粒子数N、体积V、温度T和总动量（
p,p,0）,ii为守恒量的系综（NVT）。这种情况就如同一个系统置于热浴之中，此时系统的能量可能有涨落，但系
统温度则已经保持恒定。在正则系综的MD模拟中施加的约束与微正则系综中的不一样。正则系综MD
方法是在运动方程组上加上动能恒定（即温度恒定）的约束，而不是像做正则系综的MD模拟中对运动方程加上能量恒定的约束。在正则系综MD的平衡化过程中，速度标度因子一般选下面的形式较为
合适
（4.11）
我们可将正则系综MD的Verlet算法的速度形式的模拟具体步骤列在下面：
（l）给定初始空间位置(1){t}:(i,1,2,？,N) i
(1) （2）给定初始速度{v} i
(n，1)(n)(n)(n)2 （3）利用公式r,r，hv，Fh/2m 计算在第n+1步时所有粒子所处的空间位置iiii
(n，1){r} i
8
(n，1)(n)(n，1)(n) （4）计算在第n＋1步时所有粒于的速度{v,v，h(F，F)/2m}，动能和速度标度因iiii子
（5）计算将速度n，1{v}标度因于β 值，并让该值作为下一次计算时，第n+1步粒子的速度i
n，1n，1{v,},{v}： ii
（6）返回到步骤（3），开始第n＋2步的模拟计算。
按照上面的步骤，对时间进行一步步的循环。待系统达到平衡后，则退出循环。这就是正则系综的
MD模拟过程。
下面我们举一个微正则系综的MD模拟的应用示例来看看模拟的结果。
例： 对一个总能量确定的单原于（氩）粒子系统的MD模拟计算。
我们具体选取256个原子的模拟。粒于间的相互作用位势为Lennard-Jones势
（4.12）
其中-ε为位势的极小值（取ε为能量单位），其位置在r=21/6ζ处、该体系的粒于限制在一个立方体的箱
21/2子，边界上采用最小象力约定。我们采用自然单位制。长度和时间的标度单位分别为ζ和(mζ/48ε)
-12（对氩原于该时间单位为3，10秒），这样就使得运动方程为无
量纲形式。模拟时我们考虑两个相图上的点(T*，ρ*)=(2.53，0.636)，（0.722，0.83134），它们分别具有两种立方体的尺寸，即L＝7.83和L＝6.75。初始条件假定为：各个原子处于一个面心立方格子的格点
上，而速度按相应温度下的波尔兹曼分布抽样取值。位势的截断取两个值r=2.5和r=3.6，用以比较其cc对模拟结果的影响。在执行平衡化过程中，调节粒子速度的标度因于为
(4.13)
反复上面的速度调节，直到系统能量达到给定值、在这个例子中，平衡化过程用了1000步MD模拟。模拟结果列于表4.1中。表中的误差为标准误差。系统总动能的模拟演化过程由图4.1给出。实际上，图中显示出在数百步后动能就达到平衡了。图4.2则显示出位能的平衡化过程。系统总能量的平衡化过
程则由图4.3表示，其平衡化是通过对粒子速度的调节跳跃式地达到的。图4.4为动能的分布图，模拟得到的平均速度为*＝0.3654，而理论上该值应当是。这个结果已经是相当vv,1.13T/24,0.3668不错了，因为我们只对256个粒子的系统进行了模拟。而且速度大于平均速度的粒子数所占百分比与期
望值46.7％也一致。表中的数据表明模拟结果与所选择的截断距离值变化并不灵敏。
9
表4.1 对256个粒子的氩原子系统进行1000步微正则系综MD模拟的结果
趋衡到T*=2.53，ρ*=0.636
趋衡到T*=0.722，ρ*=0.83134
图 4.1 动能演化过程图（T*＝2.53） 图 4.2 位能演化过程图（T*＝2.53）
图 4.1 总能量演化过程图（T*＝2.53） 图 4.2 动能的分布图（T*＝2.53）
10
酶动力学方法
酶学性质研究
1. 化合物对酶可逆性抑制测定
1.1 大体积稀释法
将侧活体系所需浓度的酶与10倍IC50浓度的化合物孵育30Min,以使化合物与酶充分作用。同样浓度的酶与同样体积的化合物的溶剂孵育，作为对照；孵育结束后，将上述混合物100倍稀释至反应溶液（此时，酶浓度至测活浓度，化合物浓度至0.1倍IC50)；加入底物启动反应，检测酶促-时间曲线；比较实验样品与对照稀释后的酶促曲线。
1.2透析法
将酶与高浓度的化合物孵育足够的时间，以使酶活性被充分抑制。孵育结束后，将上述混合液到大体积反应溶液充分透析，对照样品以同样的方法透析；透析完全后，测定蛋白浓度和酶活性；比较实验样品以及对照样品透稀前后的酶比活性；根据透稀前后的酶比活性的恢复情况，确定抑制剂的可逆性。
2. 化合物抑制类型研究
2.1 Km值和kcat测定
测定不同底物浓度时，酶促反应的初速度。对测得的反应速度与相应底物浓度进行回归分析，求得Km值。计算公式为米氏方程（Michaelis-Menten equation）：ν=Vmax*[S]/( Km+[S])（公式1），ν代表反应初速度，Vmax代表最大反应速度，[S]代表底物浓度。此运算可同时得到Vmax。根据公式Vmax＝kcat*[E]（公式 2）计算kcat；其中，Vmax为最大反应速度，在计算Km值时求得；[E]为酶浓度
3. 可逆性抑制剂抑制类型测定
3.1竞争型抑制剂
通过不同浓度抑制对酶的Km，Vmax 变化趋势的影响，可以初步确定抑制类型。随抑制剂浓度的增加Km增大，Vmax 基本不变。
最后，利用双倒数作图确认抑制剂性质，即在不同抑制浓度条件下，底物浓度倒数（1/[S]）对相应反应速度倒数(1/[V])做图。不同浓度的抑制剂对应的拟合直线交于Y轴同一点。
3.2非竞争型抑制剂
通过不同浓度抑制对酶的Km，Vmax 变化趋势的影响，可以初步确定抑制类型。随抑制剂浓度的增加Km基本不变，Vmax降低。
最后，利用双倒数作图确认抑制剂性质，即在不同抑制浓度条件下，底物浓度倒数（1/[S]）对相应反应速度倒数(1/[V])做图。不同浓度的抑制剂对应的拟合直线交于X轴同一点
3.3反竞争型抑制剂
通过不同浓度抑制对酶的Km，Vmax 变化趋势的影响，可以初步确定抑制类型。随抑制剂浓度的增加Km降低，Vmax降低。
最后，利用双倒数作图确认抑制剂性质，即在不同抑制浓度条件下，底物浓度倒数（1/[S]）对相应反应速度倒数(1/[V])做图。不同浓度的抑制剂对应的拟合直线平行
分子动力学方法
计算材料学
The Computational Materials Science “A new way of thinking …”
上海交通大学材料学院：金朝晖 2010.02
原子层次模拟：分子动力学方法（I）
MD simulation Observing and measuring atoms in a virtual lab
Molecular Dynamics（MD）
对一个包含N个粒子（原子或分子）的体系，给定粒子 之间相互作用势，初始条件和边界条件，通过对牛顿运 动方程做数值积分，从而得到粒子运动轨迹的方法。
五个要素： 1. 2. 3. 4. 5. 粒子间相互作用势 （interatomic potential） 初始条件（initial condition） 边界条件（boundary condition） 求解牛顿运动方程（二阶常微分方程组）的积分算法 （integrator） 粒子运动轨迹（atomic/particle trajectories）
体系的内能
对体积为V的N粒子体系，定义内能：
E ≡ K +U
其中，K为系统动能 U为系统势能
1 ? K ≡ ∑ mi x i ( t ) i =1 2
N
2
U = U ( x3N (t ))
x3N(t) 表示t时刻3D坐标集 x1(t) , x2(t) , x3(t) ,……, xN(t)
绝热体系（微正则系综, NVE）： 对于孤立体系，E为恒量，不随时间和体系所处的状态而变化 恒温体系（正则系综，NVT）：粒子数、体积、温度不变
原子间相互作用势 （1） pair potential
Lennard-Jones potential
?? σ ?12 ? σ ?6 ? V ( r ) 4ε ?? ? ? ? ? ? ?r? ? ? ?? r ? ?
Repulsion: short range Attraction : long-range Ar：σ=3.405 ?，ε = 119.8kB = 0.01032 eV (rare-gas solid: fcc packing) Electronic configuration: 1s2 2s2 2p6 3s2 3p6 (close-shell) central force interactions (no directional bonding) : good approximation for spherical atoms
初始构型：面心立方晶体（Argon）
IC (初始条件) N = 4x(5x5x5) = 500 T=0K
[001]
[010]
[100]
3维周期性边界条件 （3D PBCs）
a0=5.77?
Example
热平衡状态下的晶体
模型体系 初始条件 边界条件 原子间相互作用势 系综（积分算法） 计算目的 Ar N=500 ， fcc晶格，cubic box （Supercell） 温度 T=10 K , T=80 K 周期性边界条件 （PBCs） Lennard-Jones potential 正则系综（NVT）：粒子数、体积、温度不变 使体系达到平衡状态 观察晶体内部原子的热振动（热涨落） 以及晶体熔化过程（固—液相变）
MD Simulation
建立体系模型
t=0，给定粒子数目N，粒子间相互作用势，初始条 件（坐标、体积、温度、压力）和边界条件
平衡或弛豫
使体系在预期的温度和压力条件下达到热力学平衡
统计平均
进行足够长时间模拟，计算相关物理量
结果分析和评估
例如，径向分布函数、弹性常数
T = 10 K
T = 80 K
熔化
Hot crystal
Liquid
力的计算
rij = x ij ( t ) x= x j ? xi ij
xj(t)
? ij ≡ x
原子i 所受的力
xij rij
（单位向量）
? ?V ( r ) ? ?V ( rij ) ? ij ?∑ = fi = ?? ?x ∑ ? ? ?x i ?r r = rij ? j ≠i j ≠i ? ? 1 ?V ( r ) ? = ∑
? ? ?x ij ? r ?r r = r ? j ≠i ij ? ? (8)
原子i 和 j 之间的力 因此
fi = ∑ fij
j ≠i
? 1 ?V ( r ) ? fij ≡ ? ? ?x ? r ?r r = r ? ij ij ? ?
(9)
以及 fij = ?f ji
势能截断 truncation (cutoff)
?? σ ?12 ? σ ?6 ? = V ( r ) 4ε ?? ? ? ? ? ? ?r? ? ? ?? r ? ?
cutoff distance rc
Argon
From long-range interactions to short-range
Lennard-Jones势
rc ≈ 2.373σ
势能截断方案
?V (r ) ? V (rc ), V0 (r ) = ? ?0, r
（在r = rc，一阶导数不连续！）
?V ( r ) ? V ( rc ) ? V ′( rc )( r ? rc ), r
r
12 6 12 6 ? ? ? ? σ ?18 ? σ ?12 ? ? r ?6 ?σ ? ?σ ? ? σ σ ? ? ? ? ? ?4ε ?? ? ? ? ? + ? 2 ? ? ? ? ? ? × ? ? ? 3 ? ? + 2 ? ? ? , W (r ) = ? ?? r ? ? r ? ? ? rc ? ? rc ? ? ? σ ? ? rc ? ? rc ? ? ? ? ? ? ? ?0, ?
r
近邻原子 在 r ≤ rc 范围内的原子配位数
rc
xj(t)
1 ei = 2
rji
∑ W (r )
ji
对fcc Ar，根据左图，估算 ei和晶格常数a的关系，得到 平衡晶格常数a0，并验证此 时力（fi）是否为零。 在此基础上，如何计算晶 体的弹性模量？
the minimum image criterion: among all possible images of a particle j, select the closest (if r>rc), and throw away all the others (if r>rc) Ly
周期性边界条件 1
Lx, Ly, Ly > 2 rc
j i
rij
rc
Lx
i 的坐标 (xi, yi, zi), j 的坐标 (xj, yj, zj)
原子j的镜像
δ x= x j ? xi ij δ y= y j ? yi ij δ z= z j ? zi ij
rij =
(δ x ) + (δ y ) + (δ z )
2 2 ij ij ij
2
Maintaining neighbor lists in MD simulations
? δx ?l ? ? ij x δ xij = ? ?δ xij + l x ? ?
1 if δ xij > l x 2 1 if δ xij
周期性边界条件 2
对参照系H0，Supercell 的三个边 由向量h1， h2， h3标定，它们组 成3×3矩阵H 原子i 的坐标：
j’ i j H H0
x i = si1h1 + si 2h 2 + si 3h3 = si H 0 ≤ siμ
原子i 和 j 之间距离为：
x ji ≡ x j ? x i
如图，给定截断距离rc，在周期边界 条件下，如何判断i 和 j 是否互为近 邻原子—如果互为近邻原子，如何计 算它们之间的距离？
Supercell: one can always use 3D PBCs in MD simulations
PBCs vacuum region (>rc)
infinite system
infinite film/slab/wire Equivalent to 1D or 2D PBCs, e.g.
finite system (no PBC)
x i + lh1 + mh 2 + nh3 (l , m, n = ?∞, ∞ )
x i + lh1 + mh 2 (l , m = ?∞, ∞ )
Order-N MD Simulation With Short-Range Potential
积分算法（integrator）
t0+ L ?t t = t0
积分算法用很短的时间步长（ ?t ） 产生粒子长时间的运动轨迹。
x 3 N (t0 ) →
x 3 N (t0 + ?t ) → x 3 N (t0 + 2?t ) → ? → x 3 N (t0 + L?t )
(1) Verlet Algorithm
?t ~ 10-15 s, L~106, 粒子运动时间 ~ 1 ns
?i (t0 )( ?t ) 2 + Ο(( ?t ) 4 ) x i (t0 + ?t ) + x i (t0 = ? ?t ) 2x i (t0 ) + ? x
加速度
? f (t ) ? x i (t0 + ?t ) = ?x i (t0 ? ?t ) + 2x i (t0 ) + ? i 0 ? ( ?t ) 2 + Ο(( ?t ) 4 ) ? mi ?
?? x i (t0 ) = fi (t0 ) mi
（Taylor expansion with high order terms）
? i (t0 ) v i (t0 )= ≡x
1 [x i (t0 + ?t ) ? x i (t0 ? ?t )] + Ο(( ?t ) 2 ) 2 ?t
(2) Velocity Verlet Algorithm
To advance one step from t0 to t0+?t , x3N(t0+?t) can be obtained via
x i (t0 += ?t ) x i (t0 ) + v i (t0 )?t +
fi (t0 ) ( ?t ) 2 2mi
evaluate f3N(t0+?t), and then calculate v3N(t0+?t)
1 f ( t ) f ( t + ?t ) ]?t v i (t0 += ?t ) v i (t0 ) + [ i 0 + i 0 2 mi mi
such that one can have x3N(t0+?t) and v3N(t0+?t) simultaneously
(3) Predictor-corrector Algorithm (4) …
MD integrators
Algorithm error (intrinsic and round-off errors) stability symplectic?
… Velocity Verlet Algorithm and its high-order extensions High order predictorcorrector integrators …
The most important property for Hamiltonian systems is the Poincaré and Liuville's conservation law of phase areas…
For details see “BASIC MOLECULAR DYNAMICS” by Ju Li better long-term performance better performance at large integration time step lower computational cost with desired accuracy
Energy minimization approaches
1) Transition state : from one state to another 2) different ``attraction basins'' in the energy landscape 3) local minimum and global minimum 4) overcome energy barriers 5) relaxed state (at the bottom of the basin in the energy landscape) 6) finding optimal structure: the optimal structure would then be the one with the lowest energy Gradient descent Steepest Descent (SD) and Conjugate Gradient (CG) Methods
Approaching local minimum
? i >0 ?if fi ? x ? ? else
? i (t ) = ?i x x ? i (t ) = 0 x
原子间相互作用势 （2）
beyond L-J：嵌入原子方法(EAM)
对某个原子i，考虑它和所有其它原子之间的双体相互作用，亦考虑所 有其它原子和它的相互作用（many-body effects） 势能 = 双体相互作用（pair-wise） + 嵌入原子能量（many-body）
N 1 N U (x ) Fi ( ρ i ) = ∑V ji ( rji ) + ∑ 2 j ≠i 1 i= 3N
嵌入原子i的能量由函数F给出
ρ i = ∑ f j ( rij )
j ≠i
N
—— background density for the embedded atom i
F, f 和V 依赖具体的物理模型
(Finnis & Sinclair ,based on the tight-binding approx. Fi ( ρi ) =
ρi )
EAM一般对金属适用，尤其适合处理金属中位错、界面等问题
EAM potentials for Cu and Al
Truncation scheme of pair potential for EAM Cu
Al: valence electrons with s-p mixing Cu: valence electrons with s-d mixing EAM potentials developed by Ercolessi, Mishin et al
EAM potential: force calculations
N 1 N = U (x ) Fi ( ρ i ) ∑V
ji ( rji ) + ∑ 2 j ≠i i= 1 3N
?U ( x i ) fi = ? ∑ ?x i j ≠i
fi=
∑{
N j ≠i
rij ? ? ′ ′ ′ ′ V ( rji ) + F ( ρi ) + F ( ρ j ) ρ ( rij ) ? ? rij
}
原子间相互作用势 （3）
beyond L-J： Stillinger-Weber势
共价键 原子以共价键结合时（如Si, sp3-bonding） ，原子间的相互作 用不仅取决于原子间的距离（键长）, 而且取决于原子间的成键 方向（键角）。
势能 = 双体相互作用（two-body） + 三体相互作用（ there-body ） k
= U
i
∑ U 2 ( rij ) +
i
∑
U 3 (ri , r j , rk )
hi ( rij , rik ,θ jik ) + h j ( rjk , rji ,θ kji ) + hk ( rki , rkj ,θikj )
U 3 (ri , r j , rk ) =
i j
Tersoff- Brenner 势 （ 碳 s-p bonding）：
U
∑∑ (V
i i
R
( rij ) ? BijVA ( rij ) )
其中VR是排斥项, VA是吸引项, Bij是一个与键角θijk有关的函数
总结: MD基本概念
原子间相互作用势 （模型势的种类及特点） 周期性边界条件及应用 势函数的空间截断方案 （与原子近邻列表、周期性边界条件的关联） 积分算法 （与统计系综的关联） 热力学统计系综
ftp site of CMS lectures (updating) ftp://cms:cms@202.120.55.106
A molecular dynamics primer with examples in Fortran 90 http://www.fisica.uniud.it/~ercolessi/md/md/
LAMMPS Molecular Dynamics Simulator http://lammps.sandia.gov/
厌气混合培养中产甲烷菌和硫酸盐还原菌的动力学竞争__动力学推定的模型及实验方法
()自然科学版 中山大学学报 Vol142 No16 第 42 卷 第 6 期
2003 年 11 月 ACTA SCIENTIARUM NATURALIUM UNIVERSITATIS SUNYATSENI Nov1 2003
厌气混合培养中产甲烷菌和硫酸盐
还原菌的动力学竞争 Ξ?. 动力学推定的模型及实验方法
贾晓珊 , 李顺义
()中山大学环境科学与工程学院 , 广东 广州 510275
() () 摘 要 : 产甲烷菌 MPB和硫酸盐还原菌 SRB的竞争是一个有关厌气处理过程生死存亡的重大问题 。研
究的目的是提供一个明确实际厌气过程中 MPB 和 SRB 竞争机理的动力学研究方法 。首先 , 以 Monod 模型为基准
建立了一整套动力学模型用于微生物菌体浓度和诸动力学常数推定 。其次 , 进行包括乙酸和氢气利用 MPB 和
SRB 的两系列连续培养及批量实验 , 为下一步的动力学常数推定及讨论实际厌气过程中 MPB 和 SRB 的动力学竞
争机理打好基础 。
() () 关键词 : 产甲烷菌 MBP; 硫酸盐还原菌 SRB; 微生物菌体浓度 ; 厌气混合培养 ; 动力学模型 ; 动力学常数
() 中图分类号 : X172 文献标识码 : A 文章编号 : 052926579 20030620103204
MPB 和 SRB 之间的竞争受到各种因素的影响 。 1 前 言 一个重要的影响因素是动力学因素 。在纯培养条件
厌气过程被利用于处理生活及工业废水已有超 下获得的几乎所有的动力学常数表明 , 乙酸和氢气 过百年历史 。废水中的有机物被降解产生甲烷 —一 利用 SRB 任 何 一 种 都 能 完 全 竞 争 过 相 对 的 [ 1 ] [ 5 ,6 ] 种能量物资是此过程最大的特性。通常 , 有机物 MPB。但是 , 这些研究都是在硫酸根离子浓度 不的厌气降解是一个复杂过程 。一些聚合有机物 , 如 受限制条件下进行的 , 在硫酸根离子浓度受限制 的蛋白质 , 脂肪和多聚糖 , 由酵素首先水解形成氨基 情形下 , SRB 则不能发挥其最大的竞争能力 , 因 其[ 5 ] 酸 、脂肪酸和糖 。这些水解产物经酸化细菌形成 受到一个称为硫酸根离子饱和常数的影响 。硫 酸
() VFA 挥发性脂肪酸, 由乙酸生成细菌进一步生成 根离子饱和常数反映了 SRB 对于硫酸根离子的 亲乙酸 、二 氧 化 碳 和 氢 气 。最 后 , 由 产 甲 烷 菌 合性 , 其不但影响不同种类 SRB 之间的对硫酸 根() () () MPB经以下反应 1和 2转换乙酸和氢气 成离子的竞争 , 同时还进一步决定 MPB 和 SRB 之 间[ 1 ,2 ] 为甲烷。但是 , 乙酸和氢气也同时能作为电 子的对底物乙酸和氢气的竞争结果 。然而 , 纯培养 条
( ) ( ) 受体被硫酸盐还原菌 SRB经以下反应 3和 件下有关 SRB 的硫酸根离子饱和常数的报道非 常[ 3 ,4 ] () 4利用掉。因此 , 实际的厌气过程中 , 在MPB 有限 。
和 SRB 之间存在着对底物乙酸和氢气的竞争 。 实另外 , 与纯培养条件相比现实的条件具有很大 际上 , 如果 SRB 能有效竞争过 MPB , MPB 可能从 反的复杂性 。首先 , 这种硫酸根离子浓度不受限制的 应器完全地被排除而无甲烷生成 。这对厌气处理 过条件在实际上是很难满足的 。现实的厌气废水处理 程是一个有关生死存亡的重大问题 。 过程中 , 除了一些特殊工业污水之外 , 市政生活污 - + [ 7 ] ()CHCOO + HCH + CO 1 3 水经常只含有 50,200 mgΠL 的硫酸根离子, 因其 0’ 4 2 ΔG= - 3518 kJ 远远低于有机物浓度从而造成硫酸根离子浓度受限
()4H+ COCH+ 2HO 2 2 2 4 2 制的情形 。其次 , 现实的厌气废水处理过程中 , 绝
0’ ΔG= - 13018 kJ 非像纯培养研究那样仅仅存在单一种类的 MPB 或 - 2 - CHCOO + SO + 2H ()3 4+ - + 2CO+ 2HO 3 SRB , 而是呈现多种类 MPB 和 SRB 的混合相存 。 2 2 HS0’ ΔG= - 5713 kJ 在此 , 竞争不仅存在于 MPB 和 SRB 之间为电子来
+ - 2 - 4H+ H + SOHS源如乙酸和氢气 , 同时为硫酸根离子 , 在乙酸和氢 ()2 4+ 4HO 4 2 0’ Δ气利用 SRB 之间也存有竞争关系 。但是 , 有关在实 G= - 15212 kJ
Ξ 收稿日期 : 2003 - 03 - 28
基金项目 : 香港 RGC 资助项目
() 作者简介 : 贾晓珊 1962 年生, 男 , 教授 , 博士生导师 ; E- mail : eesjxs @zsu1edu1cn
104 () 中山大学学报 自然科学版第 42 卷 K 际的废水处理过中这些复杂的竞争关系的研究报道 s ) ( X+ P - PYΠY0 0 xp + ? μt = ln m X S+ XΠY0 0 0x很少 , 其竞争机理还远远没有被清楚地了解 。
) ) ( ( X+ P - PYΠY S + P - PY 0 0 xp0 0 p本研究的目的是提供一个明确实际的厌气废水 ln + ln X S 0 0 处理过程中 MPB 和 SRB 竞争机理的动力学研究方 ()13 法 。首先 , 以 Monod 模型为基准建立模型用于包括
( ) () () 当S+ XΠYμ K等式 12和 13分别可被简 0 0x s 微生物菌体浓度和诸动力学常数的推定 ; 进行包括
() () 化成以下等式 14和 15。 乙酸和氢气利用 MPB 和 SRB 的两系列连续培养及
) ( X+ S - S Y批量实验 。其次 , 根据建立的动力学模型及批量试 0 0 x ()μt = ln 14 m X 0验结果 , 进行包括乙酸及氢气利用在内的 MPB 和 ) ( X+ P - PYΠYSRB 的动力学常数的推定 ; 讨论实际的厌气废水处 0 0 xp ()μ15 t = ln m X 0理过程中 MPB 和 SRB 的竞争机理 。
() () 等式 14和 15为本研究的基本动力学推
2 动力学推定模型的建立 定模型 。这些模型说明 , 如果初期的底物或产物的
浓度及产率系数已知的话 , 到达某一时间的底物或 211 基本模型
产物的浓度将仅仅取决于初期的细胞浓度 。 一般 , 批量降解的生长动力学可被公式化如
下 : μ212 X和 的推定模型 0 m
d X 以两任意已知时间 t和 t的底物及产物浓度 值1 2 μ ( ) ()= - kX 5 d () () d t 分别代入等式 14或 15并两式相比 , 可得
1 到下面的等式 : S μ()d X 6 = - Y d t x) ( X+ S- SY 0 0 1 xln μX d P d S t0m 1 ()= - Y 7 p= ()16 d t d t ( ) X+ S- SY μ0 0 2 xtm 2 ln X 0μS m ()8 μ = K+ S s ) ( X+ P- PYΠY 0 1 0 xpln μμ X 这里 是比增长率 , 是最大比增长率 , K是细 μ0 m d tm 1 = ()17 ( ) μ X+ P- PYΠY t0 2 0 xp胞死亡率 , K是饱和常数 , Y是细胞产率系数 , m 2 s x ln X 0 Y是产物产率系数 , S 、 P 和 X 分别是介质中底 p
() () μ16和 17为 X和 的推定模型 。代前述已 等式 0 m 物 、产物和细胞浓度 。 () () 知的 Y和 Y入等式 16或 17则 X能被求出来 。 x p 0 () 当细胞增长在指数阶段时 , 方程 5中细胞 ()使用求得的 X代入等式 16 ()0 或 17 式的分子或分 死亡率 K是可忽略不计的 。细胞死亡率只在细胞 d μ母则 也能被求出来 。 m 增长非常缓慢和比增长率接近零时才显示出重要 [ 8 ] 213 K的推定模型 s 性。因此 , 当细胞增长在指数阶段 , 且细胞死亡
μ本研究中 , 在求出 X和 后 , 饱和常数 K 0 m s() 率 K可忽略不计时 , 如下的等式 9可从方程 d [ 3 ] () () 同样也能根据常用的 Lineweaver-Burk 图推定出 5- 8导出 :
() () 来 。等式 9和 10联立及线性化得到等式 d X d S 1 1 d P 1 - = = = μS m () X 9 - 1 Y xd t Yd t Yd t YK+ S 1 d S YKx p x s x s 1 + = ( μ( ) ) + YS- S d t S μX 0 x 0 m m当 Y和 Y两个都是恒定时 , 能获得以下等式 : x p
()18 ( ) ()X = X+ YS- S 10 0 x 0
() ()等式 18为 K的推定模型 。这里 , 以等式 18 s ( ) ()P = P+ YS- S 11 0 p 0
的左边和 1ΠS 为坐标能作成直线 , 通过直线的斜 ( ) 这里 , X, S和 P分别是初期 t = 0时介质中 细0 0 0
率 、截距及已知的常数将能惟一的求出 K。 s 胞 、底物和产物的浓度 。
() 随后 , 替代 X 或 P 入等式 9和各自进行综 3 实验方法
合化 , 分别能获得以下联合模型 :
11 富集培养 3) ( X+ S - S Y 0 0 x K s + μt = ln m ?本研究进行了如下两种类型的混合富集培养 。 XS+ XΠY 00 0x
( ) X+ S - S Y 0 0 x S (一种为仅仅含有乙酸和氢气利用 MPB 的培养 系 ()ln + ln 12 ) 列 1, 另一种则为混合含有乙酸和氢气利用 MPB X S 0 0
第 6 期 贾晓珊等 : 厌气混合培养中产甲烷菌和硫酸盐还原菌的动力学竞争 105
() 和 SRB 的培养 系列 2。两种类型的混合富集培 , 以便清除残留底物及其他杂质 。步骤为 : 起 清洗
( ) 初养分别以两个完全混合反应器 CSTR在 37 ?的 , 100 mL 的培养污泥混合液进行离心分离 , 分离
(条件下 , 使用葡萄糖 、营养素和硫酸根离子 系列 后的沉淀污泥用 200 mL 无机盐营养液进行清洗并 ) 2连续进行 。选择葡萄糖为培养底物是考虑到其 再次进行离心分离 , 分离后的沉淀污泥最终再悬浮
降解接近现实的废水处理过程 。 于 100 mL 无机盐营养液 。
() 对于添加乙酸的批量实验 表 1如 : A-1 、A- 两个 CSTR 在接种了被破碎的 UASB 颗粒污泥
AS-2 , 实验开始前 , 各小瓶首先用氮气 2 、AS-1 和 后以 人 工 底 物 开 始 连 续 培 养 。人 工 底 物 含 有 :
完全除去瓶内氧气并密封起来 , 用注射器注入 40 9 372 mg/ L 的葡萄糖 , 相当于 10 000 mg/ L 的 COD ;
mL 清洗过的培养污泥混合液 。添加浓缩乙酸底物 () 554 mg/ L 的硫酸根离子 系列 2; 包括重碳酸盐
[ 9 ] 及 p H 调整溶液以达到预定的瓶内乙酸浓度及 p H 在内的营养素和微量元素。最初 , 各 CSTR 的
() 值 。在添加硫酸根离子的批量实验 表 1如 : AS- ) (HRT 水力停留时间为 15 d 之后再减少为 10 d。
1 和 AS-2 中硫酸根离子也包括在添加的浓缩底物 ( 在各 HRT , 各 CSTR 的运行在达到稳定状态 即 ,
当中 。对于添加氢气的批量实验如 : H-1 、H-2 、 ) 连续稳定的底物去除和甲烷生产率后再连续稳定
HS- 1 和 HS-2 , 用注射器注入 20 mL 清洗过的培养污 运行 30 d , 用注射器从各 CSTR 收集适量泥样以准
泥混合液到已除去氧气并密封起来的小瓶内 , 再注 备下阶段的批量实验 。
入 2 mL 只包含无机盐的营养液或者同时包含硫酸 312 批量实验
( ) 根离子 HS-1 和 HS-2。然后 , 用注射器注入 1 atm 本研究使用前述富集培养后的系列 1 和 2 的污
( 及 20 ?的 H和 CO气 体 混 合 体 积 比 : 2 2 泥 , 在 70 mL 玻璃清液小瓶内进行动力学推定的批
) H?CO= 4?1。以上所有操作都在绝氧工作站内 2 2 量实验 。玻璃小瓶实验方法参照我们以前的方
10 进行 。 法。表 1 批量实验的条件 。对从 CSTR 收集的
培 养污泥 , 在玻璃小瓶实验之前用无机盐营养液
进行
表 1 添加乙酸的批量实验条件
Tab11 The conditions of batch experiments
实验编号 A- 1 A- 2 H- 1 H- 2 AS- 1 AS- 2 HS- 1 HS- 2
HRTΠd 15 10 15 10 15 10 15 10
添加基质 乙酸 乙酸 乙酸 乙酸 HΠCO HΠCO HΠCO HΠCO 22222222- 1 ()832 830 600 750 820 870 850 880 最初基质浓度Πmg?L - 1 ()- - - - 120 100 300 310 硫酸根离子浓度Πmg?L - 1 ()1520 2540 1570 2500 1645 1110 1640 1105 VSSΠmg?L
pH 7. 15 7. 20 7. 18 7. 10 7. 15 7. 08 7. 13 7. 20
所有批量实验在恒温震荡水浴中进行 。水浴温 然后 110 ?维持 113 min 。注射口及检测器的 min ,
温度均为 180 ?。乙酸的浓度利用同样型号的第 2 度 37 ?, 震荡强度为 35 mm ×125 strokesΠmin 。较
台气相色谱分析仪以 FID 和 10 m ×0153 mm 毛细管 强的震荡是为了实现污泥 、混合液及气相的充分混
柱测定 。升温程序为在 70 ?维持 4 min 然后 140 ? 合 , 包括促进 H和 CO从气相到液相的溶解 。批 2 2
维持 3 min。注射口及检测器的温度均为 200 ?。 量实验开始后 , 定时测量及分析生物产气量 ; 气相 以 40 mL/ min 流速的氦气作为载气 。硫酸根离子浓 中甲烷 、H及 CO的浓度 ; 混合液中的 p H、乙酸 2 2 () 度利用离子色谱分析仪 岛津 HPLC-10A测定 。 及硫酸根离子浓度 。
参考文献 : 313 分析方法 VSS 和 p H 的测定利用标准方法 。生物气体中 [ 1 ] MCCARTY P L , SMITH D P. Anaerobic wastewater
氢气和甲烷的百分比用装备了 TCD 检测器的气相 () treatment[J ] . Environ Sci Tech ,1986 , 20 12: 1200 -
( ) 色谱分析仪 HP-5890 ?测定 。毛细管柱为 25 m 1206.
[ 2 ] SPEECE R E. Anaerobic biotechnology for the industrial ×0153 mm 的 CarboPLOT P7 系列 。以 30 mL/ min 流
() wastewater treatment[J ] . Environ Sci Tech ,1983 ,179: 速的氩气作为载气 。升温程序为 50 ?维 持 215
106 () 中山大学学报 自然科学版第 42 卷
416 - 427. and growing conditions[J ] . Arch Microbiol ,1984 ,137 : 26 [ 3 ] WIDDET F , PFENNIG N. A new anaerobic , sporing , acetate-- 32.
oxidizing , sulfate- reducing bacterium , desulfotoma- culum [ 7 ] YODA M , KITAGAWA M , MIYAJ I Y. Long term
( ) emendacetoxidans [ J ] . Arch Microbiol , 1977 , 122 :119 competition between sulfate- reducing and methane-
- 122. producing bacteria for acetate in anaerobic biofilm[J ] . Wat [4 ] BADZIONG W , THAUER R K. Growth yields and growth Res ,1987 ,21 : 1547 - 1556.
( ) rates of desulfovibrio vulgaris Marburg growing on [8 ] YANG S T , OKOS M R. Kinetic study and mathematical
hydrogen plus sulfate and hydrogen plus thiosulfate as the modeling of methanogenesis of acetate using pure cultures of
sole energy sources[J ] . Arch Microbiol ,1978 , 117 : 209 - methanogens[J ] . Botech. Bioengrg. 1987 ,30 :661 - 667.
214. [9 ] J IA X S , FANG H H P. Methanogenic characteristics of [5 ] INGVORSEN K, ZEHNDER A J B , JORGENSEN B B. formate- utilizing sludge [J ] . J Environmental Engineering ,
() Kinetics of sulfate and acetate uptake by desulfobacter 1999 , 125 7:596 - 601.
postgatei[ J ] . Appl Environ Microbiol , 1984 , 47 : 403 - 408. [10 ] IA X S , FURUMAI H , FANG H H P. Extracellular [ 6 ] ROBINSON J A , TIEDJ E J M. Competition between sulfate- polymers of hydrogen- utilizing methanogenic and sulfate-
reduction and methanogenic bacteria for Hunder resting 2 reducing sludges[J ] . Wat Res ,1996 , 30 : 1439 - 1444.
Kinetic Competition bet ween Methanogenic and
Sulfate- reducing Bacteria in Anaerobic Mixed Cultures
?: Model Development
JIA Xiao-Shan , LI Shun- yi
()School of Environmental Science and Engineering , Sun Yat- sen University , Guangzhou 510275 , China
Abstract : Bacterial mass concentrations and kinetic parameters of acetate- and hydrogen- utilizing methane-producing
() () bacteria MPBand sulfate- reducing bacteria SRBwere estimated to examine the kinetic competition between MPB and SRB in anaerobic mixed cultures. The cultures were enriched by two chemostat reactors , one was fed with glucose alone
while the other was fed with both glucose and sulfate. The parameters were estimated from batch experiments data , based
on the Monod model and developed kinetic models.
( ) ( ) Key words : methane-producing bacteria MPB ; sulfate- reducing bacteria SRB ; bacterial mass concentration ; anaerobic mixed culture ; kinetic parameter
