分子模拟时间越长越好吗？ - 分子模拟 - Gromacs/Amber/NAMD

15楼2011-07-01 11:23:12

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子，希望有用哦~

对接的分子用于分子动力学模拟的问题已经有12人回复
分子模拟中MSD的问题讨论已经有25人回复
扩散系数的计算及体系平衡的判断已经有4人回复
分子动力学的时间已经有9人回复
怎样看gromacs模拟后的小分子和蛋白质之间的相互作用的残基和氢键？已经有14人回复
gromacs结果分析怎样分析氢键能和结合能已经有12人回复
LAMMPS中ReaxFF力场是怎么回事？已经有17人回复
蛋白质构象模拟已经有12人回复
关于材料的热力学稳定性结构稳定性已经有16人回复
金属功函数的测量方法有几种已经有5人回复
什么动力学模拟软件可以模拟化学反应过程已经有13人回复
如何评价误差和变化？已经有10人回复
做离子液体的分子模拟，用哪个分子模拟的软件好呢？已经有3人回复
gromacs 中g_msd问题已经有8人回复
有没有人用Amber做过膜蛋白动力学模拟已经有14人回复
生测数据的标准差大于平均值，怎么办？已经有7人回复
AMBER中qm/mm的问题已经有5人回复
【讨论】分子模拟还是同源建模已经有16人回复
【求助】msd算法已经有10人回复
【求助】NAMD和Gromacs哪个简单？已经有5人回复
【求助】关于用分子动力学来做光谱的问题已经有14人回复

点击这里搜索更多相关资源

科研从小木虫开始，人人为我，我为人人

wavingsea9911

木虫 (小有名气)

模拟EPI: 1
应助: 2 (幼儿园)
金币: 638.9
帖子: 231
在线: 73.7小时
虫号: 565734

★ ★ ★ ★ ★ ★
奋斗1s(金币+1):谢谢参与
zh1987hs(金币+5, 模拟EPI+1): 谢谢 2011-07-01 19:49:20
奋斗1s(金币+4): 恩，谢谢您的指点！ 2011-07-01 21:55:33
奋斗1s(金币+1): 2011-07-01 22:02:18

只对一个蛋白质构象进行的动力学模拟，一般用于对同源建模的结果进行构象优化，目标是去除模型中的高能键、高能角等不合理的构象。我假定楼主也是这个目的。那么时间是否越长越好呢？
  我个人认为不是这样。因为根本的一个问题就是，作为动力学模拟起点的模型构象未必是（或者说一定不是）该蛋白质的真实构象。以这个起点开始模拟过程，随着时间的推移，我们期望的结果是：该构象会向真实构象靠拢，最后稳定在真实构象上。但实际的情况更可能会是：各种构象偏差逐渐放大，最后稳定于一个完全不同于实际的构象，或者整个结构发生坍塌，发展成一团不知所以的肽链。
  另一个问题是溶剂体系的问题，如果在真空条件下对酶蛋白，特别是膜蛋白进行动力学优化，如果不加约束，很难得出正确的结果。
  但值得注意的是，错误的动力学模拟结果的蛋白质评价指标也可能会很好，因为动力学模拟一定会让体系向能量最低的方向发展。所以，蛋白质评价指标只是佐证，不是铁证。
  最后的结论是，用于蛋白质模型优化的动力学模拟，时间长短不是很重要，个人认为1-5ns都可以。关键是要对该蛋白的真实构象有一个概念，以此控制动力学模拟的走向，保证系统向正确的方向发展。
  8楼的阐述很细致，是很好的理论基础。遇到具体问题还要具体分析。

赞一下(5人)

17楼2011-07-01 15:31:55

ChemiAndy

木虫 (正式写手)

模拟EPI: 40
应助: 52 (初中生)
金币: 2573.3
帖子: 590
在线: 553.6小时
虫号: 836907

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+10, 模拟EPI+1): 谢谢CHEM 2011-07-02 23:28:30

引用回帖:

Originally posted by wavingsea9911 at 2011-07-01 02:31:55:
只对一个蛋白质构象进行的动力学模拟，一般用于对同源建模的结果进行构象优化，目标是去除模型中的高能键、高能角等不合理的构象。我假定楼主也是这个目的。那么时间是否越长越好呢？
我个人认为不是这样。因为 ...

感谢17楼补充。我在8楼说了很多，其实没有很好地回答楼主的问题，做Docking，模拟时间要多久。

Docking的本质是看看哪个药物小分子（底物）与蛋白的结合最好。我不知道楼主打算从哪几个方面来考察对接是否“好”，好不好有很多定性和定量的指标。最常用的热力学指标有平均能量差（焓变），和结合自由能。不过，最重要的还是要看你的底物和结合位点的功能关系。比如有的结合位点的功能是要打开蛋白的通道的，那对接后通道尺寸就是主要指标；有的位点是催化反应的，那么就要求结合构型有利于反应，结合自由能再大也不行滴。

从楼主模糊的描述来看，似乎有能量或自由能（“对接结果”）和蛋白构型（“蛋白质评价”）两种评价标准。而且后2ns不如前3ns，这里面说明2个问题，1.不同的评价指标有不同的弛豫周期，对模拟时间的要求不一样；2.前后模拟结果有变化，说明体系尚未达到平衡，或者说已经平衡，但是是在单个弛豫周期内单调变化。

建议：1. 明确衡量对接评价指标以及相应的热力学/几何结构的微观观测量；2. Plot这些微观观测量，看看是否已经平衡，是否已经完成一个弛豫周期，如果没有，就一定要继续探索下去。3. 看些相关体系（相同蛋白或同家族蛋白）Docking方面的Review，了解这个领域的通常的问题，要求和广泛接受的做法。

最后回应一下17楼wavingsea9911 的不同意见，17楼从一个特定的例子 -- 同源建模以后对体系的松弛（结构优化）为例 -- 说明模拟时间长短并非越长越好。这个我是同意的。并非只有以采样为目的的模拟才是模拟，所以模拟时间不能一概而论。

但是17楼说模拟的构象会坍塌到一个“错误”的构型（“实际的情况更可能会是：各种构象偏差逐渐放大，最后稳定于一个完全不同于实际的构象，或者整个结构发生坍塌，发展成一团不知所以的肽链。”）这个我觉的稍有不妥。相空间中没有错误构型，只有能量高的构型。既然我们说相空间，这里就包含了一切可能构象。只是有些能量太高，出现概率小到可以忽略。但是既然这部分是高能构象，那体系怎么可能“坍塌”到高能构象上呢？这有3个可能：1力场错了，2这个高能构象恰好是你初始构象的一个局部最小，与全局最小之间的能垒很高，滑到里面出不来了；3这个构象跟你说的“最稳定”的折叠构象的能量很接近，它也属于这个蛋白在水里面出现概率很大的一种状态。因此，这3种情况不能用来证明“模拟时间未必越长越好”。只要是在相空间进行的有效采样（力场正确），都给我们关于蛋白构象分布的有用信息。当然，这里我只是在字眼和大道理上较真，其实17楼只是在强调采样必须有效effective（“控制动力学模拟的走向，保证系统向正确的方向发展。”）而不是否定采样必须充分sufficient。不过鉴于这些概念对于初学者很重要，我就补充一下。

[ Last edited by ChemiAndy on 2011-7-2 at 10:10 ]

23楼2011-07-02 22:29:07

yongleli

木虫 (正式写手)

模拟EPI: 1
应助: 93 (初中生)
金币: 3787.6
帖子: 742
在线: 99.9小时
虫号: 303595

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
奋斗1s(金币+1):谢谢参与
zh1987hs(金币+10, 模拟EPI+1): 非常感谢 2011-07-05 12:47:28

引用回帖:

Originally posted by ChemiAndy at 2011-07-02 22:29:07:
感谢17楼补充。我在8楼说了很多，其实没有很好地回答楼主的问题，做Docking，模拟时间要多久。

Docking的本质是看看哪个药物小分子（底物）与蛋白的结合最好。我不知道楼主打算从哪几个方面来考察对接是 ...

这里我要回答ChemiAndy版友的一个小疑问，并讲另外一个关于模拟时间的问题。

关于傅献彩的物化，我们本科时候学习的是另一本，比较偏向化工，习题很多是学校的老师们在做项目时遇到的问题。虽然公式没怎么详细推导，但是让我们觉得学了物化就能干很多事了。我的毕业后在化工设计院工作的同学们确实也是觉得受益匪浅。后来到了南大和山大考察，发现原来老师们虽然用傅献彩的书，但是并不按照书本教。能否把具体的熵焓等概念讲明白就看老师个人功力了。一般是可以的。这也和两校物化课时特别多有关。

做分子动力学模拟，是否模拟时间越长越好？如果计算的模型能够正确反映真实的体系，我觉得是时间越长越好。毕竟理论上，时间系综经过无穷长时间可以遍历相空间中每一个点。但是实际上我们都知道这是不可能的，而且有时候初条件会严重影响模拟的结果。第二点，比如要模拟binding过程或者steered MD，一条轨迹往往也不够。因为如C版友所言，事件的发生有几率性，而MD由于模拟的时间有限（真实世界中1s是很快的时间，但是全原子分子模拟中这几乎是不可能的任务），一条轨迹未必能在有限时间中观测到体系发生想要的变化。

但是还有一点需要注意，MD模拟两大基本问题一是抽样，另一个是力场。
目前的力场我们一般认为它们是行为很好的，可以反应生物大分子的性质。但是某些情况下，比如带电配体的binding，含金属的生物分子，传统的力场很可能不能很好的模拟。用传统力场模拟的话，短时间内体系跟XRD结构接近，但是长时间之后就会变成非自然的构象。比如上世纪90年代Kollman模拟了一个体系，只做了几百ps，但是后来的研究发现这个体系用AMBER力场来做的话在4ns之后会变成非自然构象。这是因为固有的力场中能量最小点不是天然产物的能量最小点。这方面的工作比如
http://pubs.acs.org/doi/abs/10.1021/ja807374x
揭示了用传统的力场本质上是不能正确模拟的，遑论抽样了。

赞一下(2人)

26楼2011-07-05 12:37:12

zh1987hs

金虫 (著名写手)

★ ★ ★ ★ ★ ★
奋斗1s(金币+1):谢谢参与
ghcacj(金币+5, 模拟EPI+1): 谢谢，这个帖子很热啊 2011-07-05 13:35:09

非常感谢楼上诸位高手的讨论~作为一个MD的纯粹应用者来说，理论上很深刻的东西还是不懂，就我们应用的角度来说，主要考虑的是MD能为我们的实验理论解释做一些什么，就我们进行的工作来说，主要思考以下的问题：
1.热稳定性，能否通过MD找到蛋白的热敏感区域，从而对其进行针对性的改造
2.活性位点关键氨基酸的作用，通过长时间的MD,找到蛋白和配体结合时可能发生的构象改变，找到特定氨基酸对结合影响的原因
那么结合上个专家所言，时间和力场应该是结合起来考虑的，如果力场不合适，比如现在没有发现能够很好模拟金属酶的力场(QM/MM除外)，那么进行长时间的模拟也不一定能取到蛋白质的真实构象，采用现在通用的AMBER,CHARMM或GROMOS的话，不知道力场对结果的影响大不大，在这种情况下，我觉得长时间的模拟还是很有必要的，现在很多文献都做到ms级别，用于研究酶的变构效应，毕竟许多反应需要的时间不是ns级别的，可能需要将时间延长到一定尺度才有可能观察到对应的变化，一点拙见，请指正

赞一下(2人)

27楼2011-07-05 12:54:30

ChemiAndy

木虫 (正式写手)

模拟EPI: 40
应助: 52 (初中生)
金币: 2573.3
帖子: 590
在线: 553.6小时
虫号: 836907

★ ★ ★ ★ ★

送鲜花一朵
zh1987hs(金币+5, 模拟EPI+1): 谢谢 2011-07-06 08:01:50

To yongleli大侠:

1，你所给文献的例子揭示的应该是模拟前“验证力场准确性”的重要性，而不是从根本上否定非极化力场。

该文献说是经典非极化力场不能反映蛋白-配体结合位点的一根关键氢键，然而使用极化力场就行，模拟与实验一致。俺觉得这不能称为“用传统的力场（非极化）本质上是不能正确模拟”，因为力场方法和所有的实验方法一样，不同的力场有其特有的精度precision和准确度Accuracy两方面。精度是说其采样的事件的能量可以精确到多大，准确度是相对平均值的偏差RMSD(误差)。采样精度由力场决定。采样准确度由模拟方法以及事件本身波动大小（不确定性）决定。

(2)力场的精确度precision决定所能模拟什么样的体系，性质，和过程

非极化力场的精确度在10kcal/mol的数量级（5~20kcal/mol）。极化力场可达1kcal/mol数量级。因此，简单说，如果一个体系，它两种状态的能量差在10kcal/mol，我们就能用力场来计算它在这两种状态下的分布概率的差。

具体到蛋白构型，单根氢键H-bond的能量恰在～10kcal/mol，在你所给文献中，非极化力场不能正确描述它完全有可能。

一般情况下，对于模拟蛋白构型，非极化力场还是能够胜任的，因为蛋白构型折叠与非折叠构型的能量差显然比较大，解折叠的过程的能垒也比较大。已发表的很多模拟已经告诉我们，非极化力场在模拟构型方面与实验的符合度是基本令人满意的。然而，如果要比较折叠状态下，个别位点的几何结构，这就取决于该位点不同结构的能量差，如果小于10kcal/mol，那么由于力场精度问题，得出的结果可能就是错误的。

再然而，错误也是有限的。想象一下我们拿了一个尺子，如果它的精确度是10cm，意味着测量值的十位上数字是精确的，个位上估计的，个位以下没有意义。假如我们用这个尺子比较一堆木棍的长度，那么，只要所有的木棍长度大于1cm，那么尽管我们不能准确知道每根木棍的长度，但是其相对长度还是能够依靠个位数字上的估计值比较得出。这恰如自由能的计算，自由能是用力场这根尺子，去测量每个构象的能量，从而得出在每种构象上的分布概率，再由概率确定各种性质，或者过程前后的热力学自由能差和动力学自由能垒。也就是说，如果构型的能量差虽然小于10kcal/mol，但是如果大于1kcal/mol，结果仍然具有一定意义。

上面力场精确度我讲的是数量级，因为不同力场质量也不一样。怎么办呢？模拟前要验证力场准确性。这也是前贴的重点。

总之，在当前非极化力场水平下，若该性质由能量小于10kcal/mol数量级的能量波动确定，则定量预测实验值没有意义，但是，定性结果基本满意。

(3) 力场精确度对模拟时间的影响

对平衡态性质，模拟时间由该性质的弛豫周期决定；对非平衡过程的话，模拟时间由该过程的能垒决定。

[ Last edited by ChemiAndy on 2011-7-5 at 11:32 ]

赞一下(3人)

29楼2011-07-05 23:42:45

ChemiAndy

木虫 (正式写手)

模拟EPI: 40
应助: 52 (初中生)
金币: 2573.3
帖子: 590
在线: 553.6小时
虫号: 836907

★ ★ ★ ★ ★
zh1987hs(金币+5, 模拟EPI+1): 谢谢 2011-07-06 08:01:38

To zh1987hs版主：

（1）确定蛋白热稳定性的“热敏感区域”，其实就是确定unfolding过程的起始步。假定存在多个竞争的起始unfolding区域，考虑到每个起始步有不同的能垒，那么，我们可以模拟多次来确定从每个区域起始的概率。但是，如果这些起始unfolding步的能垒差较大，在10kcal/mol左右，预测结果应该没问题。能量差小了，恐怕就不好说。不过，用可极化力场可以大大提高结果准确性。

所需的模拟时间呢，由起始unfolding能垒确定，代入阿伦尼乌斯公式，可以得到特征反应时间t，取这个时间的3-10倍，就是你的模拟要跑的时间尺度。当然，如何估计起始unfolding能垒就是个大问题。一般由实验确定，或者类似体系的模拟结果推测。

（2）考察活性位点残基与底物结合的几何构型变化所需的模拟时间，也是一个能量问题，即这些构型变化索要翻越的能垒。方法同上。准确性分析也类似。

30楼2011-07-06 00:03:35

ChemiAndy

木虫 (正式写手)

模拟EPI: 40
应助: 52 (初中生)
金币: 2573.3
帖子: 590
在线: 553.6小时
虫号: 836907

★ ★ ★ ★ ★

送鲜花一朵
zh1987hs(金币+5, 模拟EPI+1): 谢谢 2011-07-06 08:02:09

引用回帖:

Originally posted by gzgzgz at 2011-07-05 13:34:06:
"非极化力场的精确度在10kcal/mol的数量级（5~20kcal/mol）。极化力场可达1kcal/mol数量级"
这个估计是如何得到的？针对某些测试体系的benchmark后的回归统计结论？是大部分的极化力场还是某种 ...

问的好！

这个还真没有出处，印象的来源是导师和读过的一些paper。硕士的老板，是开发力场的，这么提过。后来读博士选的一个分子模拟的课，那老师是Mackerell的学生，发展CHARMM和可极化力场的，他也大致这么说过，加深了印象。我感觉他们也未必能肯定的说个结论。只能说，这个估计是个数量级估计，是根据分子模拟能否适用于特定体系和过程的经验总结，而非严密的统计结论。举个实例，液体混合焓一般在正负几个kcal/mol，俺一哥们用了3-4年的时间搞清了一件事，排除模拟方法的影响，所有非极化力场不能得到任何定性或定量准确预测结果；极化力场（induced dipole model）可以得到定性准确的预测结果（不同体系和组成）。几年青春得到的结论阿，兄弟们。再比如，离子水合，有没有极化，得出的water-water相关时间差别一倍以上。TIP3P（非极化）的扩散系数是实验值的2倍。采用极化力场以后，这些性质和实验值的差别大大缩小。而对于这些性质，恰好是在零点几个到数个氢键的精度范围，所以说10kcal/mol与1kcal/mol分别是非极化力场和极化力场的精度范围，大差不差。

这里你可能会问，不同力场之间，质量差别大不大？还是经验结论，非极化力场这种粗糙的模型，模拟凝聚态体系性质，主要是vdW参数和静电作用参数起作用。只要这些参数的取值在一定的合理范围内，所得到的模拟结果差别不是很大。为什么呢？一般含极性分子或离子，金属等的体系，vdW的吸引能作用并不大，主要起一个排斥作用，决定体系的体积性质；凝聚体系之所以能凝聚，还是静电参数决定的，即体系的电荷分布，或者说介电常数。如果介电常数差别不大，其中的带电粒子的运动差别就不大。对于有些精细的过程，能量差在几个kcal/mol的，不同的力场得出的结论不一样，纯属偶然，没有任何意义。非极化力场根本不足以表现这些细节过程。

[ Last edited by ChemiAndy on 2011-7-5 at 16:42 ]

32楼2011-07-06 04:54:25

ChemiAndy

木虫 (正式写手)

模拟EPI: 40
应助: 52 (初中生)
金币: 2573.3
帖子: 590
在线: 553.6小时
虫号: 836907

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+10, 模拟EPI+1): 谢谢！ 2011-07-06 15:16:08

报道的自由能“精确到”2-3kcal/mol，应该是准确度Accuracy，即独立采样结果间的均方偏差。衡量的是多次采样值在平均值周围分布的偏离程度。这个与该模型系统的波动（不确定度）和采样方法即FEP技术有关。举例来说，比如体系的尺寸很小，那么波动就大。波动值大概与粒子数目的平方根的成反比sigma~1/SQRT(N)。另外，FEP方法中，你的Window（微扰步数）取的越多，你的结果就越“精确”。这个报道的精度值必然有其计算方法的参考文献，你可以去验证一下。

上面说的模拟的数值结果的精确度Accuracy是2-3kcal/mol，与其本身物理精度precision只有10kcal/mol并不矛盾。因为分子模拟的数值精度可能达到小数点后很多位，它们虽然不是物理意义上准确的能量值，但却对数值模拟的稳定性(能量守恒)很重要，对衡量模型体系的尺寸效应，采样的有效性很重要，我们当然可以按这些“高精度”的数值结果计算出来达到个位数的数值误差。

问题是，用一把精度为10kcal/mol的尺子，“估算”出自由能在2-3kcal/mol的采样波动，是否有物理意义？这个很少有论文去讨论这个报道的准确度与力场精度本身的联系。我的想法是，自由能衡量的是什么呢？是体系在不同能量状态上的分布概率。只要产生的概率分布是对的，自由能就可能是精确的。那精度为10kcal/mol的尺子，我们能估读能量到个位，则必然能反映由个位数能量决定的概率分布的不同，从而得出2-3kcal/mol的精确度。这个估读有物理意义，则自由能准确度也有物理意义。

因此，教科书上提到FEP都说是“估算”自由能，这个估算指的就是物理意义上的估算，就是指用低精度模型（low precision）去大致估量可能超出当前物理精度范围的过程的自由能过程。而所报道的“精确度”(Accuracy，or error, 就是“结果+正负2～3kcal/mol”) 指的是数值意义上的精确度。

Precision与Accuracy定义与区别：精度Precision是采样结果的有效位数，是由测量工具即“尺子”决定的；准确度Accuracy(or say error)是各独立采样结果相对平均值的偏差程度，与体系本身性质的波动和采样方式有关。前者属于绝对误差，后者属于相对误差。想象一个问题，具有同样数值位数的ab initio计算结果与分子力学计算结果有啥差别？

[ Last edited by ChemiAndy on 2011-7-6 at 02:13 ]

34楼2011-07-06 15:06:56