M11T6.58120.482J 生物系统中的算法与引技术果田 Rruce Tidor教授 Jacob K.White较授 MT6.581/20.482问题集2 2006.3.14,犀期二，5pm 1,建立一个抽获高子的动态模型，这里我们将探讨几种分子动力学中的积分 方法。图】显示了一个这样的系统一一个移动离子被一个方阵的类带电离子 抽获， oblle lon 忍1分了司小学的物具仿直承统 Malb脚木文件moldyn.m包含了对这个系统的描述，同时静电相互作用的 力专能量方程和van der Waals相互作川方程如下： 9qi 4mEorii (1) 2) van der Waals参数e和a表示最有利的相互作用(e)和使相互作用变得 不再有利的使离()。下表列出了3个不同离子的参数。 lon m (g/mol) q(e) G(A) e [kcal/mol) F- 18.998 -10 2.73295 072000 d- 35.453 -10 4.41724 0.11779 Br- 79.904 -1.0 4.62376 009000
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 1 MIT 6.581/20.482J 问题集 #2 2006.3.14，星期二，5pm 1. 建立一个捕获离子的动态模型。这里我们将探讨几种分子动力学中的积分 方法。图 1 显示了一个这样的系统——一个移动离子被一个方阵的类-带电离子 捕获。 图 1 分子动力学的物理仿真系统 Matlab 脚本文件 moldyn.m 包含了对这个系统的描述，同时静电相互作用的 力与能量方程和 van der Waals 相互作用方程如下： (1) (2) van der Waals 参数ε和σ表示最有利的相互作用(ε)和使相互作用变得 不再有利的距离(σ)。下表列出了 3 个不同离子的参数
MT6.58120.482J 生物系统中的算法与计算技术基础 Bruce Tidor教授 ace的K.White教授 请注意，在分子级别的仿真中，距离用M1×10),电荷单位为质子电荷 (e,质量单位是g/mol,能量单位为kcal/mol:在这些单位下，ra是2397×10e2 A-(kcal/mol)。虽然选择这样的单位是的时间单位必须为4.888821×10s,但 通常都将时间单位转化为ps(在atlab脚本文件moldyn.m中实现). ()实现前向做拉积分法。并将离子的动力以0.01ps的步长传播。 (b)描述移动离子的动态轨迹 ()系统能量（运动的，电势的和总能量）如何跟随时间变化？ (d山使用0.001s的步长重复运行仿真程序，并重复b)和（©），解释他们的不 同 (e)实现Velocity Verlet积分法，并将离子的动力以0.O1ps的步长传播. 重复b)和(c)。如何测试这个轨迹的有效性？验证你的观点。 (f)用0.1ps的步长重做(d),并解释区别. (g)采用C1和Br作为移动高离子（选择你认为最合适的积分法）重复仿真。将 轨迹与以F作为高子的仿真结果比较，离子的轨迹怎样随着离子的类型变 化？离子的能量又怎样变化？ ()初始位置为9.0A的情况下再次运行仿真，同样变换各种离子类型。讨论 轨迹和能量的观察结果。 2.远程静电交互作用。在这个问题中，我们将使用特征值分解。通过TLAB的 eg函数，说明有时简单的近似能够显著的降低计算复杂度（空间和时间），同 时不需要牺牲太多的准确性。对一个矩阵A,其特征值入和相应的特征向量x满 足： Ax (3) 而对称矩阵有正交的特征向量完备集。 图2显示了两个点电荷线。每个原子上都有部分电荷，我们将其表示成位于 原子中心的点电荷。若用向量可表示簇1中的原子电荷集：5表示簇2中的原子 电荷集。g中电荷在簇2中产生一个电势场：用。表示簇2中原子中心的电势向 量，由？产生的电势向量可由矩阵乘法得到 =Pq (4) 其中P是电势矩阵：P的第1列是簇1中位置1的单位电荷对能2中原子中 心产生的电势集。（如此线1中的任意电荷分布？分别对每个电荷产生一定量的 回应。)我们已经建立了一个对称的物理系统，这样P就是对称的。一般情况下
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 2 请注意，在分子级别的仿真中，距离用 Å(1×10-10m)，电荷单位为质子电荷 (e)，质量单位是 g/mol，能量单位为 kcal/mol；在这些单位下，ε0 是 2.397×10-4e 2 ⋅Å⋅(kcal/mol) -1。虽然选择这样的单位是的时间单位必须为 4.888821×10-14 s，但 通常都将时间单位转化为 ps（在 matlab 脚本文件 moldyn.m 中实现）。 (a) 实现前向欧拉积分法，并将离子的动力以 0.01ps 的步长传播。 (b) 描述移动离子的动态轨迹 (c) 系统能量（运动的，电势的和总能量）如何跟随时间变化？ (d) 使用 0.001ps 的步长重复运行仿真程序，并重复(b)和(c)，解释他们的不 同。 (e)实现Velocity Verlet 积分法，并将离子的动力以0.01ps的步长传播。 重复(b)和(c)。如何测试这个轨迹的有效性？验证你的观点。 (f)用0.1ps的步长重做(d)，并解释区别。 (g)采用Cl- 和Br- 作为移动离子（选择你认为最合适的积分法）重复仿真。将 轨迹与以F- 作为离子的仿真结果比较。离子的轨迹怎样随着离子的类型变 化？离子的能量又怎样变化？ (h)初始位置为-9.0Å的情况下再次运行仿真，同样变换各种离子类型。讨论 轨迹和能量的观察结果。 2. 远程静电交互作用。在这个问题中，我们将使用特征值分解，通过 MATLAB 的 eig 函数，说明有时简单的近似能够显著的降低计算复杂度（空间和时间），同 时不需要牺牲太多的准确性。对一个矩阵 A，其特征值 λ 和相应的特征向量 x 满 足： Ax=λx (3) 而对称矩阵有正交的特征向量完备集。 图 2 显示了两个点电荷簇。每个原子上都有部分电荷，我们将其表示成位于 原子中心的点电荷。若用向量 q 表示簇 1 中的原子电荷集；s 表示簇 2 中的原子 电荷集。q 中电荷在簇 2 中产生一个电势场；用ϕ 表示簇 2 中原子中心的电势向 量。由 q 产生的电势向量可由矩阵乘法得到 ϕ = Pq (4) 其中 P 是电势矩阵；P 的第 i 列是簇 1 中位置 i 的单位电荷对簇 2 中原子中 心产生的电势集。（如此簇 1 中的任意电荷分布 q 分别对每个电荷产生一定量的 回应。）我们已经建立了一个对称的物理系统，这样 P 就是对称的。一般情况下
MT6.58120.482J 生物系统中的算法与计算技术基础 Bruce Tidor教授 acbK.White教授 P不是对称的，因比通常采用奇异值分解(SVD)来进行分析而不是特征值分解. Trefethen和Bau的书数字线性代数中对SVD有详细的说明. 镜透 电而镇 图2静电系统：再个分散的，对称的电待线 如果簇很大，炬阵P也相应的变大，其存储和做乘法运算的开销就会变大。 请下载脚本文件escm并用你喜欢的任何编辑器打开，这个脚本中有许多空白 行，你可自己加入详细代码： (©)填写矩阵P:元素（化，应该等于位于源点j的单位电荷所产生的到目标点 1的静电电势。方程(1)已经给出了1和)之间相互作用的能量，该能量等 于1的电荷乘以)到1的电势。检查P是否是对称的。 (b)使用MATLAB的©g函数将P分解为特征向量矩阵V和特征值对角矩降 D,检查 P-VDVI (5) 并且V的列是正交的（即VV,=1当-少，或0其他）.因此-（换 句话说是霄矩阵)。当分开的距离很小和很大时，观察D中元素，其 范围（特征值集合）如何变化？ (©)假设D中某一个对角线元素正好为0。则Pg在相应特征向量对应的方向 上没有值，如果你愿意可以自己验证一下：推荐到Trefethen和Bau的书 中了解详细信息，同样的，如果一个特征值非常小，比如说是0Fa) (计算机的精度)，则P在（对应的特征向量）这个方向的组成值就 很可能非常小 如果我们能够提前知道，我们就能使用一些主特征值，所有值都会变得 很小，可以采用以下方法来设计对P的近似：当特征值足够小时去掉所 有的特征值特征向量对。接着，用低阶近似P米乘，则Pg变成 Pq=Pq=VDv'q (6)
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 3 P 不是对称的，因此通常采用奇异值分解（SVD）来进行分析而不是特征值分解。 Trefethen 和 Bau 的书数字线性代数中对 SVD 有详细的说明。 图 2 静电系统：两个分散的，对称的电荷簇 如果簇很大，矩阵 P 也相应的变大，其存储和做乘法运算的开销就会变大。 请下载脚本文件 estatic.m 并用你喜欢的任何编辑器打开。这个脚本中有许多空白 行，你可自己加入详细代码。 (a) 填写矩阵 P：元素(i, j)应该等于位于源点 j 的单位电荷所产生的到目标点 i 的静电电势。方程(1)已经给出了 i 和 j 之间相互作用的能量，该能量等 于 i 的电荷乘以 j 到 i 的电势。检查 P 是否是对称的。 (b) 使用 MATLAB 的 eig 函数将 P 分解为特征向量矩阵 V 和特征值对角矩阵 D，检查 P=VDV-1 (5) 并且 V 的列是正交的（即 = 1 j T Vi V 当 i=j，或 0 其他）。因此 V-1=VT （换 句话说 V 是酉矩阵）。当分开的距离很小和很大时，观察 D 中元素，其 范围（特征值集合）如何变化？ (c) 假设 D 中某一个对角线元素正好为 0。则 Pq 在相应特征向量对应的方向 上没有值。如果你愿意可以自己验证一下；推荐到 Trefethen 和 Bau 的书 中了解详细信息。同样的，如果一个特征值 λi非常小，比如说是 O(εmachine) （计算机的精度），则 Pq 在 vi（对应的特征向量）这个方向的组成值就 很可能非常小。 如果我们能够提前知道，我们就能使用一些主特征值，所有值都会变得 很小，可以采用以下方法来设计对 P 的近似：当特征值足够小时去掉所 有的特征值/特征向量对。接着，用低阶近似 Pˆ 来乘，则 Pq 变成 Pq Pq VDV q T ≈ ˆ = ˆ ˆ ˆ (6)
MT6.58120.482U 生物系统中的算法与计算技术基健 Bruce Tior教授 aec动K.White教授 其中广表示对应于k个最大特征值的特征向量，D是这些特征值组成的 ×k的矩阵。当={1,2.8}计算P。 (@改变簇之间的分码距离，绘制距离增加时的相对误差Pg-户 ()假设每个族中有：个电荷。存储P需要多少内存空间？需要多少次浮点 操作才能计算得到P? ()如果不存储P而存储其k阶的近似P,需要多少内存空间？需要多少次 浮点操作才能计算得到户g?如果用趋近于一个非常大的值，P和P之间 的比较又会怎样？
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 4 其中Vˆ 表示对应于 k 个最大特征值的特征向量， Dˆ 是这些特征值组成的 k×k 的矩阵。当 k={1, 2, 8}计算 Pˆ 。 (d) 改变簇之间的分隔距离，绘制距离增加时的相对误差 Pq Pq − ˆ (e) 假设每个簇中有 n 个电荷。存储 P 需要多少内存空间？需要多少次浮点 操作才能计算得到 Pq？ (f) 如果不存储 P 而存储其 k 阶的近似 Pˆ ，需要多少内存空间？需要多少次 浮点操作才能计算得到 Pˆ q？如果 n 趋近于一个非常大的值，P 和 Pˆ 之间 的比较又会怎样？