计算化学公社

标题: 有关两相体系模拟的问题请教 [打印本页]

作者
Author: kunkun 时间: 2015-1-19 15:56
标题: 有关两相体系模拟的问题请教
sob老师你好，
由于课题需要模拟两相体系的蛋白构象变化，但是在构建过程中产生一些问题，想请教一下！
使用的程序是gromacs，蛋白力场是AMBER SB99 ,有机小分子使用的是GAFF
根据JUSTIN的两相体系教程，http://www.bevanlab.biochem.vt.e ... sic/01_genconf.html。
构建好box，并且insert-molecules，-nmol设定为1000，填充过量后（此处的过量是指无法继续填充）实际填充807分子。。跑em和NVE平衡。结果跑散了，出现了真空。
我自己分析是由于随机插入分子，并没有达到实际饱和的状态。
(, 下载次数 Times of downloads: 129)

经过思考，有几个问题实在困扰。

1 命令gmx insert-molecules和gmx genconf，两个填充命令我都尝试过，insert会跑出真空，并且有向两层结构发展的趋势。如果我预先填充好有机小分子，
再经过一定时间的平衡，应该也会出现两层体系，此时产生的gro文件也无法直接利用。不知justin的教程中，为何没有跑散。
而genconf在em时，结构会散架，产生游离的氢和碳链小片段。（请问这个是否正常？是否是由于pbc条件产生的）。
sob老师是否有更好的建模方案推荐？使得从一开始建模就是有机相比较饱和的形态？packmol初步了解输入文件是pdb，会不会不好应用GAFF和AMBER力场？
附上EM的mdp (, 下载次数 Times of downloads: 3)

2 根据justin的教程，是预先把有机相体系跑了NVE和10ns的NPT平衡后，得到的坐标文件，才进行填充水分子，
但是过程中把范德华半径从0.1改到0.3~0.4 ，根据我的尝试，需要改成0.7才无水分子填入有机相，但这样改，会对范德华力的计算产生误差。
是否需要把半径在mdrun前，修改回去？

3 使用acpype产生的itp和gro文件，我拷贝了力场ff ，nonbonded.itp中加入了新原子类型。在后续的run中，会不会对结果有比较大的影响。
应为GAFF力场atomtype是小写。（其实根本问题是，GAFF和AMBER力场在gromacs中的应用是否合适，因为我看有人发过相关的文章）

作者
Author: sobereva 时间: 2015-1-19 17:00
1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

用packmol建模完全可以，比gmx自带的强大得多。这和用什么力场完全无关。

你也可以先单独模拟很小的纯的有机分子盒子，跑平衡了，也就很致密了，加有机分子的时候直接把这个有机分子盒子往真空区域填充。（加水时候用的spc216.gro就是前人预先平衡好的216个水分子的盒子）。相反，也可以把有机盒子跑平衡后把盒子尺寸扩展，再往空着的区域加水。

2 这和计算过程无关，但还是应当及时改回去，免得以后搞乱。
不改范德华半径也无妨，可以直接用VMD载入文件，在保存新结构的时候，用范围选择语句把指定部分的水给去掉，这样更灵活准确。

3 只要操作过程正确就完全适合。有大写有小写这是必然的，别搞乱就行。

作者
Author: kunkun 时间: 2015-1-19 20:38

sobereva 发表于 2015-1-19 17:00
1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

瞬间就明白了，感谢sob老师。缩小体系跑平衡再扩充效果好很多！也帮助了理解gromacs的gro概念，
我继续演研究packmol 对比建立模型的质量哪个好。再次感谢！

作者
Author: kunkun 时间: 2015-1-24 19:44

sobereva 发表于 2015-1-19 17:00
1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

sob老师您好，我目前做动力学用的是gromacs 但是在一些设置的问题上有点疑问，希望老师指点一二。我经常看到有文献记载Electro- static interactions was modelled using a coulombic cut-off set at 8 A ̊ and 14 A ̊ .
我想请问一下这里的coulomb cut-off为何有2个值？但是rcoulombic只能设定为一个值.此处是否指的是rcoulombic和rvwd的值呢？是否和开关函数和PME有关系？我阅读gro手册时它推荐的却是rc=rw=rlist..

作者
Author: sobereva 时间: 2015-1-24 21:24

kunkun 发表于 2015-1-24 19:44
sob老师您好，我目前做动力学用的是gromacs 但是在一些设置的问题上有点疑问，希望老师指点一二。我经常 ...

应该是twin-range cutoff，即8埃以内的静电相互作用每一步都重新计算，而8~14埃的静电相互作用是每隔nstlist步才重新更新一次。
对于cutoff方式计算范德华作用，rlist和rvdw对应这两个值。
对于cutoff方式计算静电作用，rlist和rcoulomb对应这两个值。

若用PME方式计算静电作用，rlist必须等于rcoulomb，没法用这种twin-range

作者
Author: kunkun 时间: 2015-1-24 21:38
本帖最后由 kunkun 于 2015-1-24 21:43 编辑

sobereva 发表于 2015-1-24 21:24
应该是twin-range cutoff，即8埃以内的静电相互作用每一步都重新计算，而8~14埃的静电相互作用是每隔nstl ...

sob老师，我对您解释的这个问题还是有疑问，和我之前的理解不太相同，按照我的理解rlist应该是临近列表的原子对，计算短程相互作用的力（二面角，键伸缩等)，不是一般设置为0.9~1nm么?
若设置rc=rw=rlist时，所有的力都在此计算，在此外的静电力按照PME实空间和虚空间的加和，
当rc和rw大于rlist是双截断，在rc和rw的力都是每一步计算，而在rlistlong在rc和rw之间的力是NS期间计算一次。rlistlong之外的静电力按PME计算。
不知是否严重理解错误...
按照老师所说，那在rcoulombic 之外的静电力不就产生了误差了么（即使配合开关函数）。。我还一直理解为PME是配合双程截断使用的....现在对rlistlong的理解，是否作用就在于配合开关函数，在rlistlong处力矫正为0？？还是rlistlong的存在是为了缓冲区间？

作者
Author: sobereva 时间: 2015-1-24 21:54

kunkun 发表于 2015-1-24 21:38
sob老师，我对您解释的这个问题还是有疑问，和我之前的理解不太相同，按照我的理解rlist应该是临近列表的 ...

你的理解有误。

rlist、cutoff这些都是对范德华、静电这两类非键相互作用而言的，bond,angle,dihedral,1-4那些和成键有关的作用与此没有丝毫关系。

当库仑和范德华作用都使用cutoff时，rlist小于rcoulomb和rvdw时是twin-range，即rlist内的作用每步都计算，rlist与rcoulomb和rvdw之间的是每nstlist步更新一次。更远部分的作用完全忽略。

PME计算静电作用的情况与此截然不同。只有一个cutoff，rlist必须等于rcoulomb，在rcoulomb以内的在实空间计算，在rcoulomb以外的在倒易空间计算。都是每一步都计算，而且这样处理静电作用是精确的，不管多远的静电相互作用都会被计算，因此PME是计算周期性体系静电相互作用最理想的方法。而且rcoulomb在理论上只影响计算效率而不影响计算结果。

switching function在一般的研究中没用。PME时静电作用本来就是精确的，对于范德华作用，由于收敛速度远比静电作用快，rvdw设得大一些也就基本计算精确了。

作者
Author: kunkun 时间: 2015-1-24 22:09

sobereva 发表于 2015-1-24 21:54
你的理解有误。

rlist、cutoff这些都是对范德华、静电这两类非键相互作用而言的，bond,angle,dihedral ...

原来如此!上次我还给别人解释错了，马上去纠正...感谢sob老师指导

作者
Author: kunkun 时间: 2015-1-24 22:17
本帖最后由 kunkun 于 2015-1-24 22:20 编辑

sobereva 发表于 2015-1-24 21:54
你的理解有误。

rlist、cutoff这些都是对范德华、静电这两类非键相互作用而言的，bond,angle,dihedral ...

The MD simulations were conducted with NAMD (28) and performed with periodic boundary conditions at a constant pressure (1 atm) and temperature (300K) using Langevin dynamics with Nosé-Hoover Langevin piston pressure control. The long-range electrostatic interactions were calculated using the particle mesh Ewald algorithm. Both the electrostatic and the Lennard-Jones interactions had a twin-range cutoff of 10-12 A.
老师，那这个文献内德，为何静电力PME可以和twin-range cutoff 一起使用.........难道是软件之间对此概念的设计不同么？

哎呀，懂了，老师当我没说。。哈哈

作者
Author: kunkun 时间: 2015-1-24 22:57

sobereva 发表于 2015-1-24 21:54
你的理解有误。

rlist、cutoff这些都是对范德华、静电这两类非键相互作用而言的，bond,angle,dihedral ...

sob老师我对老师说的还有个小疑问，是否可以存在这种设置方法？
rlist=rcoulombic 使用PME计算，而在计算范德华力时，可存在twin-range？我看有些文献就是这么设置的，不知是否正确？

作者
Author: sobereva 时间: 2015-1-24 23:16

kunkun 发表于 2015-1-24 22:57
sob老师我对老师说的还有个小疑问，是否可以存在这种设置方法？
rlist=rcoulombic 使用PME计算，而在计 ...

这样没问题。我比较习惯rlist=rcoulomb=rvdw=1.2

作者
Author: kunkun 时间: 2015-1-25 00:27

sobereva 发表于 2015-1-24 23:16
这样没问题。我比较习惯rlist=rcoulomb=rvdw=1.2

谢谢老师我明天试试用这个跑

作者
Author: kunkun 时间: 2015-1-27 19:41

sobereva 发表于 2015-1-24 21:24
应该是twin-range cutoff，即8埃以内的静电相互作用每一步都重新计算，而8~14埃的静电相互作用是每隔nstl ...

sob老师你好，我在做动力学的时候结果很好，我想提取出某一帧的pdb文件后继做分子对接。
但是当时为了节省计算时间，把box调整地刚好。但是在跑的过程中蛋白分子移动了（即使我用了移除重心的平动。）正好穿过了pbc。我在vmd中提取pdb的时候，分子被砍成了两半= =
我又尝试了gmx traconv ，-pbc中选择whole，但是处理完后就无法和gro坐标对应上了...直接读取trr也没有图像显示...
请问sob老师有没有更好地办法解决这个？
（我看手册中pbc有一个periodic-molecules的选项，但是他的解释是分子会分开输出.我就没敢用）

作者
Author: kunkun 时间: 2015-1-27 19:42
本帖最后由 kunkun 于 2015-1-27 21:55 编辑

sobereva 发表于 2015-1-24 21:24
应该是twin-range cutoff，即8埃以内的静电相互作用每一步都重新计算，而8~14埃的静电相互作用是每隔nstl ...

找到解决办法了..

作者
Author: kunkun 时间: 2015-2-8 13:05

sobereva 发表于 2015-1-19 17:00
1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

sob老师，我在两相体系的模拟中出现了几个小问题，想向老师请教一二。
我模拟的是以十六碳烷和十碳烷近似长链和中链的脂肪酸对酶构象变化影响，但是在模拟中十六碳的影响不合符实验结果（构象变化幅度太小，多次重复的结果也很不稳定，有时候有变化有时又无变化）
我仔细查阅了一些资料，我怀疑引起的原因是由于两相界面体系压力控制选择的问题，引起界面压力
我参照了一篇膜蛋白的教程作为我的体系近似。

1 发现其中各向异性的环境模拟，压力控制使用的是semiisotropic或anisotropic，但我之前的模拟使用的是等方性的控压，而有机相的压缩比较大。（不知对密度也有了大的影响，在gmx energy中无法定义新的组单独检测有机相密度或压力..)，现在对界面压力的选择方式比较迷茫....
ps：之前模拟结果不稳定，我还以为是由于小分子电荷优化不足引起..调试多次计算后，发现影响并不是很大。

2 在膜蛋白的教程中，使用了模拟退火缓慢升温和升压，不知道我所模拟的体系是否也需要模拟退火，更好地平衡界面与蛋白之间的接触？稳定界面的体系。

3 我所研究的方向文献的模拟方法阐述的太简单，我对一个小的细节一直不清楚如何处理。
蛋白质的模拟，初始的构象是晶体，但这只是蛋白的一个瞬间，我们是否有必要在界面模拟之前，先将蛋白在水中模拟足够长的时间，待其稳定后再随机选择其中某一个构象，重新能量最小化，在进行界面体系的建模？

问题比较多..望老师见谅

作者
Author: sobereva 时间: 2015-2-8 16:18
我不清楚你模拟的体系具体是什么样的，最好有个截图。

2 压力一直保持1atm就行了。缓慢升温可以做，这是比较谨慎的做法。

3 这个我想没太大必要。

作者
Author: kunkun 时间: 2015-2-8 16:47

sobereva 发表于 2015-2-8 16:18
我不清楚你模拟的体系具体是什么样的，最好有个截图。

2 压力一直保持1atm就行了。缓慢升温可以做，这是 ...

我的体系是以以下截图为初始模型，请sob老师帮我看看
(, 下载次数 Times of downloads: 56)

界面是沿着xy轴变化的。
每次模拟时长是10ns秒，一般文献上构象变化在2ns左右，我做了大概20多个重复。只有2次是发生了构象变化。
不知是否是由于界面不够稳定和蛋白作用的点，有关系？

我仔细观察过使用十碳链的有机分子，对蛋白构象作用比长链的强（收敛的比较快~1ns）
但是实验的结果是长链的效果更强..

作者
Author: sobereva 时间: 2015-2-8 20:29
蛋白质的朝向是怎么确定的？不同朝向下模拟结果可能有较大不同。

构象发生变化，是相对于什么状态而言变化的？水环境下的？如果是这样的话，先在水下模拟可以试试。

此体系用semiisotropic控压应该比较理想。至于你说的界面压力问题，这我不好说。

也注意考虑控温的影响。不知你现在是各个相单独控温还是怎么做的，可以尝试改改控温设定看看如何产生影响。有时把温度提高点可以加速构象变化，更早、更容易地出现预期的现象。

作者
Author: kunkun 时间: 2015-2-8 21:05
本帖最后由 kunkun 于 2015-2-8 21:24 编辑

sobereva 发表于 2015-2-8 20:29
蛋白质的朝向是怎么确定的？不同朝向下模拟结果可能有较大不同。

构象发生变化，是相对于什么状态而言变 ...

蛋白质温度因子波动比较大的几条螺旋链侵泡在有机相中，使得酶转变为活性状态，通过gmx editconf中rotate选项调整方向，方向我一直都是固定的。水环境下的模拟已经做过了，rmsd 0.2nm 没有构象变化..

压力控制semiisotropic，1bar、tau_p 1.0 xy压缩率4.5e-5 z 压缩率为0 （之前z压缩也是4.5e-5，结果体系沿着z轴过度压缩，周期性边界的蛋白距离小于2nm...、现在尝试z轴压缩为0.但是这情况下压力收敛比较慢）不知此处设置是否适合。
并且，comm-grps = protein_c16 （有机和蛋白作为一个移除组） sol_NA

我现在的控温是V-rescale，阴阳离子、蛋白、有机相、水，都是单独控温的，温度308k tau_t 0.5 。不知老师为何提出单独控温这个问题？（老师是指单独调高有机相体系的温度？）

作者
Author: sobereva 时间: 2015-2-8 21:24
设定倒是没什么问题。

主要是怕你没有单独控温。
可以所有组的温度都提升比如20K看看有没有什么新情况

作者
Author: kunkun 时间: 2015-2-8 21:28

sobereva 发表于 2015-2-8 21:24
设定倒是没什么问题。

主要是怕你没有单独控温。

好的！谢谢sob老师，我先去尝试尝试！

作者
Author: kunkun 时间: 2015-2-9 13:19
本帖最后由 kunkun 于 2015-2-9 13:25 编辑

sobereva 发表于 2015-2-8 21:24
设定倒是没什么问题。

主要是怕你没有单独控温。

sob老师，我在使用 semiisotropic 空压进行NPT时，发现如下两种设置，对模拟体系产生的结果都不是很满意。

pcoupl             = Parrinello-Rahman
pcoupltype       = semiisotropic
tau_p             = 5.0
ref_p             = 1.0       1.0
compressibility = 4.5e-5       4.5e-5
这一种设置，z轴单独控压时，进行10ns模拟中，盒子被不断拉伸，最后蛋白脱离有机层...
但这种设置压力在平衡时比较快达到预期值。（1bar）

pcoupl             = Parrinello-Rahman
pcoupltype       = semiisotropic
tau_p             = 5.0
ref_p             = 1.0       1.0
compressibility = 4.5e-5       0
而这个设置，xy轴并没有拉伸，z轴被我固定。但是在npt的平衡中（3ns），压力
一直在60bar周围波动（我设定的时1bar）...而盒子的形状没有多大的变化，（我预计x y轴应该是会拉伸才对，毕竟z轴被固定了）

这个控压方式貌似也不太适合我的蛋白体系。难道我还是用回isotropic?

还是附上我的mdp吧

(, 下载次数 Times of downloads: 7)

(, 下载次数 Times of downloads: 5)

(, 下载次数 Times of downloads: 1)
我根据教程make_ndx 蛋白和有机相为一个组，SOL和阴阳离子一个组，移除重心平动。
是否这个设置只是适用于膜蛋白而非两相有机体系的蛋白模拟呢？（毕竟蛋白是固定在模中，而我的体系是需要蛋白在有机层中发生激活的..比较大的构象变化）

作者
Author: sobereva 时间: 2015-2-9 15:02

kunkun 发表于 2015-2-9 13:19
sob老师，我在使用 semiisotropic 空压进行NPT时，发现如下两种设置，对模拟体系产生的结果都不是很满意 ...

控压肯定是得semiisotropic。
我建议先模拟个纯粹的水+有机相体系，不含蛋白，等这个体系模拟平衡了，盒子尺寸也不变了，压力都稳定到1par了，再把这个盒子进行适当地扩展（如果需要的话），把蛋白塞进去模拟。
如果直接带着蛋白模拟，盒子尚未稳定，一旦盒子出现拉伸蛋白也容易跑出去。

作者
Author: kunkun 时间: 2015-2-9 15:12
本帖最后由 kunkun 于 2015-2-9 15:14 编辑

sobereva 发表于 2015-2-9 15:02
控压肯定是得semiisotropic。
我建议先模拟个纯粹的水+有机相体系，不含蛋白，等这个体系模拟平衡了， ...

我昨晚睡前也想过这个办法，但是不知道怎么才能把蛋白填入一个饱和的盒子中。用vmd的位置选择也只能是x y z的横向操作..而且我的蛋白一半是侵泡有机相当中.

作者
Author: sobereva 时间: 2015-2-9 15:22

kunkun 发表于 2015-2-9 15:12
我昨晚睡前也想过这个办法，但是不知道怎么才能把蛋白填入一个饱和的盒子中。用vmd的位置选择也只能是x y ...

VMD选择功能极为灵活，挖个球很容易，比如
within 5 of xxx就代表选择xxx附近5埃内的，以外的就是exwithin
还可以比如sqr(x-5)+sqr(y+4)+sqr(z) > sqr(5) 选择在(5,-4,0)的半径5埃以外的

作者
Author: kunkun 时间: 2015-2-9 15:28

sobereva 发表于 2015-2-9 15:22
VMD选择功能极为灵活，挖个球很容易，比如
within 5 of xxx就代表选择xxx附近5埃内的，以外的就是exwith ...

谢谢sob老师，我这就去建模试试！

作者
Author: kunkun 时间: 2015-2-9 16:11
本帖最后由 kunkun 于 2015-2-9 16:33 编辑

sobereva 发表于 2015-2-9 15:22
VMD选择功能极为灵活，挖个球很容易，比如
within 5 of xxx就代表选择xxx附近5埃内的，以外的就是exwith ...

sob老师，我刚才操作了一下，感觉挖个球感觉这个做法不太可行，因为我的蛋白分子有500多个氨基酸。挖个球基本都占掉了整个体系的一大块了（球不够大的话，用gmx inser-molecules，塞不进去分子）。而且这样操作对蛋白质的方位也不好确定。
我个人一个的思路：
我之前是单纯使用有机相进行了10ns的NPT平衡。（使用的时gmx insert-molecules 设定个数为1000 box 5 5 5 ，可能填充的过量了）
但是由于蛋白体系中的水出现，导致有机相和水的模拟和真空中模拟的情况不同，实际的密度不一样，
导致两相体系box的z轴发生急剧的膨胀。
于是我想先把有机相和水在一定条件下模拟平衡后，再用vmd除去水分子，调整有机相的盒子范围。
再建立一个调整好方向的蛋白体系。使用gmx solvate 把处理后的有机相填充。调整范德华半径，再填充水。
此时再去平衡，密度上的温度应该就解决了？（在经过能量最小化，NVE npt的平衡，以此去近似已经平衡好的两相体系）

不知老师觉得方案是否可行？或则其他的建议？

作者
Author: sobereva 时间: 2015-2-9 16:44
蛋白很容易塞进去。最好的办法是，先模拟水+有机相，平衡后，如果不够大就往XY两方向各延展一次。把蛋白的pdb里的坐标标复制到此体系的pdb文件里。在VMD里用mouse-move-molecule平移、旋转蛋白到合适的位置。然后在VMD里保存成新的pdb文件时选all not same resid as (within 1.5 of protein)就行了，和蛋白重叠的分子就都去掉了（这和把蛋白塞进膜里操作一样，更多说明见http://hi.baidu.com/sobereva/item/e1c7874c603580adde2a9fda）

你的做法也可以一试，但最理想而且最简单的还是上面这种。

作者
Author: kunkun 时间: 2015-2-9 16:52

sobereva 发表于 2015-2-9 16:44
蛋白很容易塞进去。最好的办法是，先模拟水+有机相，平衡后，如果不够大就往XY两方向各延展一次。把蛋白的p ...

谢谢sob老师，我先研读一下您的教程！

作者
Author: kunkun 时间: 2015-2-10 15:33
本帖最后由 kunkun 于 2015-2-10 16:57 编辑

sobereva 发表于 2015-2-9 16:44
蛋白很容易塞进去。最好的办法是，先模拟水+有机相，平衡后，如果不够大就往XY两方向各延展一次。把蛋白的p ...

sob老师，我按照您的教程思路，并根据实际情况，改善了一下我的做法。但出现了一些错误

我的做法：因为z轴被压缩过大（从7nm增至21nm），一部分的水和有机相通过vmd位置选择去除后另存为pdb（已经same residue）。然后genconf 把已经经过1ns平衡（达标1bar）的模型沿着xy轴扩增1次（-nbox 2 2 1），打开vmd另存为pdb。由于是两相体系扩增，水和有机小分子在pdb中的序号是交替穿插的，于是我手动调整，他们在pdb中的行号，然后gmx editconf重新编辑原子序号，ultraedit再手动调整分子序号，再修改top文件。完成了模型的矫正。grompp也没有任何问题。

后来mdrun时出现了，过多的LINCS WARNING..然后我发现他们的坐标和盒子的大小不对称。于是我用vmd，通过7 move 调整了坐标和盒子。但是警告依然出现。（虽然可以设置环境变量，但是我担心结果不对，所以没做）
再后来我发现这种方法的能量过高了，尽管能量最小化后的能量，Potential Energy = 3.2940355e+20...如此高德能量体系应该崩溃....
后来我想着按照box的压缩率去设计初始的盒子大小和分子数，但是这样做并不是很稳定也太随机，模型多的话，探索也比较耗时..就舍弃了这个方案。请教一下sob老师这个是我哪部出了错？

我认为是在扩充盒子那一步应该出错了，不知老师所说的xy轴扩充一次的命令是否是genconf?还有在vmd中调整坐标和box大小的关系...还比较迷糊。望老师解答一下。手动调整貌似太草率了一点，我认为这步也可能出了点问题..

————————————————————————————
#后来发现这种方法需要editconf 重新定义box 和center ，然后再按照sob老师的膜蛋白构建教程，把蛋白放入体系，重新能量最小化，再走一遍NVE NPT 。。。那之前的平衡岂不是白费...

作者
Author: sobereva 时间: 2015-2-10 17:50
没搞懂你的意思。

Z轴自发增至21nm后，应该各方向都已经平均处在1bar了。参考这个变化过程，你用相同的方法建立一个初始体系，恰当选取盒子尺寸和填充区域，使得它经过模拟恰好最后差不多Z的尺寸合适，并且水和有机分子占一半一半（如果一次做不到这么巧，可以多试几次）。然后用genconf把它们在X、Y方向平移复制（如果原先盒子尺寸不够而需要这么做的话）。然后，再模拟一下新得到的这个体系体系，确认一下是否确实盒子尺寸不变了，而且压力也维持好了。如果没问题了，那么这个水+有机分子的盒子就很理想了。之后再把蛋白插进去。

你之前弄得模拟体系大抵是因为原子间有明显不合理接触导致崩溃。

作者
Author: kunkun 时间: 2015-2-10 18:01

sobereva 发表于 2015-2-10 17:50
没搞懂你的意思。

Z轴自发增至21nm后，应该各方向都已经平均处在1bar了。参考这个变化过程，你用相同的 ...

明白老师的意思了。
主要是z轴我控制得比较小，大概在7~8nm之间，按照原来的压缩率，z轴起始得设为2..
留给有机小分子的距离大概就只有1nm....水分子1nm....

我之前的意思是，模拟z轴偏大一点（比如到15nm），然后用vmd除掉上面和下面多出来的体积，然后以这个体系进行模拟。结果崩溃了

作者
Author: sobereva 时间: 2015-2-10 19:44
1nm也无妨。
反复试，总能找到合适的设定。

还是不清楚你的意思。但怀疑你出问题的原因是与相邻盒子相接处出现了原子不合理的接触。最好在VMD里把镜像也显示出来看看是否有这种情况。

作者
Author: kunkun 时间: 2015-2-10 20:04

sobereva 发表于 2015-2-10 19:44
1nm也无妨。
反复试，总能找到合适的设定。

哈哈哈谢谢sob老师，找到问题了。
应该是有原子的碰撞了，我在之前的基础上把体系x y z 用vmd裁剪掉之后，扩充xy轴，直接加入蛋白，重新能量最小化，再走一遍2次平衡，盒子也基本不变。观察压力，均值范围在4~-4 左右。我看差不多就直接MD了。结果果然比之前的好太多了。不到100ps就开始构变了

作者
Author: kunkun 时间: 2015-2-13 03:14
本帖最后由 kunkun 于 2015-2-13 03:31 编辑

sobereva 发表于 2015-2-10 19:44
1nm也无妨。
反复试，总能找到合适的设定。

1 sob老师，我在阅读文献的时候，发现别人模拟的有机小分子比如辛烷等。最终经过一定时间的md后，会聚集收拢，形成一个微粒状的球体。
但是我许多次的模拟中发现，有机小分子最终会形成一个圆柱形的长条（由于pbc存在，分子在收拢过程中一旦穿过边界就会连接到另一头的分子，最后变成圆柱...)
我怀疑过是不是由于我体系的水太少了，距离不够，但是我对照文献的方法，按照他们的建模（对应他们的水数量和有机分子量）。最后还是形成圆柱条了.....

比如这是我阅读的一片文献中，他们最后md的图：
(, 下载次数 Times of downloads: 97) (, 下载次数 Times of downloads: 77)

2  请问下sob老师，这种体系的控压方式仍然是选择semiisotropic？（体系的起始状态仍然是有机层和水层），我觉得对于这种interface的模拟还是isotropic更优...

3  还有个问题。是关于这种体系下，comm移除重心组的设置应该是：蛋白有机分子 SOL_NA 、还是蛋白_有机分子  SOL_NA 、亦或者是蛋白 non-protein、system....
对于这个概念不是很清楚...

4  还有..求老师证实一下现象有无问题.. 在NVE平衡的时候，设置-DPOSE 。另外是否需要同时限制水分子的位置呢（加上-DPOSE_WATER)
若除了蛋白以外都不加上限制，蛋白和水在NVE的过程中，会发生收缩（有机小分子的收缩更加明显）产生一定的真空。（不知道这个现象是不是正常的?毕竟NVE时压力都已经快-600bar）

5 关于能量最小化mdp设置时，对于蛋白质，按照老师的经验，emtol的设置一般是多少为好呢？（我看国外的论坛，有人说蛋白设置一般是1000即可，若做自由能计算的话，需要更小。）

这几天埋头苦做苦看文献..把之前很多的小疑问都一并打上来了..

感觉做蛋白构象变化，都是蛋白质的某一个瞬间...能不能成功，感觉随机性挺大,,,观察多次md轨迹有感。

作者
Author: sobereva 时间: 2015-2-13 11:27
1 这和PBC有关。如果随意模拟一堆辛烷比较随意地分散在水中，那么最后肯定是变成一个球。但如果你模拟一层辛烷的膜，膜又比较厚，那么这个膜会一直维持。如果膜不是很厚，恰好水分子在一侧穿越了膜，和膜另一侧的水联通上，那么膜就会收缩成柱。如果水分子的运动又把这个柱切断了，那么柱就会收缩成球。

2 如果已经收缩成球了，就用isotropic。如果是柱，应该semiisotrpic，让柱的方向控压和另外两方向不同。

3 蛋白_有机分子 SOL_NA 可能合适一些

4 如果都做了限制就没什么意义了。其实没必要在NVE下做平衡。毕竟初始结构中烷烃很难堆得很紧密，不控压的话跑起来出现真空是正常的。

5 这个不重要，只要MD能跑起来就行了。很多情况，觉得初始构型比较好，我都不做优化而直接跑（但需要经历0K升温过程，直接耦合到实际模拟温度容易崩）

MD的随机性太大了，前一阵子模拟个膜，结果很好，和实验很相符，结果重新算一次来验证，结果却相反了，立刻心情低落。

作者
Author: kunkun 时间: 2015-2-13 14:47

sobereva 发表于 2015-2-13 11:27
1 这和PBC有关。如果随意模拟一堆辛烷比较随意地分散在水中，那么最后肯定是变成一个球。但如果你模拟一层 ...

谢谢sob老师！
同感啊，心情和结果挂钩..做了30多个模型，跑出了2个..也算欣慰了

最后一天在实验室，希望做出点结果好回家过年！

作者
Author: kunkun 时间: 2015-2-16 02:35
本帖最后由 kunkun 于 2015-2-16 02:40 编辑

sobereva 发表于 2015-2-13 11:27
1 这和PBC有关。如果随意模拟一堆辛烷比较随意地分散在水中，那么最后肯定是变成一个球。但如果你模拟一层 ...

sob老师，我还想请问一下关于模拟退火升温的两个小问题...

1 我在模拟中设置的时single，起始温度是273K 升温至 313K ，time是0~300ps 后来在模拟中我发现，体系的起始温度是263，比我设定的起始温度还要低10度。再升温过程中，的的确确是刚好升温加热40度。但是还没达到目标水平..很是奇怪，是否是因为我起始的温度和速度没有在NVE下产生好的缘故？然后我续跑AN。结果一开始体系就直接飙到313K了..

2 在模拟退火的过程中，是否需要设置速度的产生呢？我看国外论坛答案不一，有人提出的方案是现在NVE下模拟50ps后产生初始速度，再AN缓慢加热（NPT），进行压力和温度的同时平衡。也有人直接AN就设置温度的产生和直接NPT。老师认为哪个方案比较好呢？AN缓慢升温的同时进行控压，是否已经可以取代传统NVE升温+NPT控压的平衡方案了？（生物体系）
ps：
3 我在模拟酶在有机相的构象变化时，发现对于有机相的控压，预先isotropic进行控压1ns，待体系的各向压力都趋于稳定时，MD时在采取semiisotrpic控压，box的大小会控制的比较好。老师认为我这种自创的方案是否可行？会不会对体系有比较大的偏差？

4 我在多次的MD模拟中，发现swiss同源建模（同源99%）出来的模型，然后用save服务检测，模型质量都比modeller要好...但md的结果却很随机（失败居多），后来我使用数据库的蛋白结晶pdb，clean之后，md，反而结果重复性高很多（基本都能重现）。于是我猜想这是由于建模的问题引起的，老师对与同源建模swiss评价如何？（我尝试过modeller的单模板和多模板，感觉可控性都不是很强，saves检测模型质量也没有swiss的评分高，不知老师有无其他的推荐？
还是说这种情况我直接就用pymol进行点突变？（我大概要联合同时突变3~5个氨基酸）

5 关于gromacs GPU加速的问题，我想请问下老师，我在网上查询到openMM库可以使用GPU加速，但是我在编译的时候只有GPU选项而没有openMM，我是否需要下载GPU版的gromacs？那这个适合我的单机么？同时我也不太清楚openMPI和这两个有什么关系..希望老师能解答一下（我们实验室只给我配备了一台i7+gtx650的单机..不是专业计算，资金也有限...）我想能不能有什么办法可以提高我的计算效率！

预祝各位sob老师，新年快乐！

作者
Author: sobereva 时间: 2015-2-16 14:22
1 grompp之后会显示升温从多少到多少，看看是否设错了。这和之前跑NVE没关系。如果实在搞不定，就在正式模拟前再跑一小段先平衡到313K。

2 没必要做NVE，直接升温就行了，初始速度随机。

3 没什么问题，反正最终的状态和设定都是无误的。

4 你都有了结晶pdb为何还同源模建？如果你要模拟的和已有的pdb只是几个残基不同，直接突变就行了，用不着同源模建。

5 以前的gmx才需要openMM，而且限制很多，现在的gmx自己就直接支持GPU加速，用不着openMM。gmx不分CPU和GPU版，只取决于编译时是否带着GPU的选项。openMM和openMPI没有丝毫关系。GTX650太low了，弄个GTX770不到两千块钱就已经能比单纯用CPU跑快几倍了，是穷人跑MD的首选。这里有一些相关讨论
单路计算化学攒机配置推荐
http://hi.baidu.com/sobereva/item/4960f63ac5ab3c99b90c0329

作者
Author: kunkun 时间: 2015-2-16 15:43
本帖最后由 kunkun 于 2015-2-16 15:44 编辑

sobereva 发表于 2015-2-16 14:22
1 grompp之后会显示升温从多少到多少，看看是否设错了。这和之前跑NVE没关系。如果实在搞不定，就在正式模 ...

谢谢sob老师，由于我的蛋白结晶里面，有非标准的残基修饰和糖，所以我想着同源建模更方便直接....
sob老师我还想问一个问题就是如果我使用gromacsA54a7 糖分子和非标准残基应该我自己用高斯用HF算还是有别的在线服务器可以生成？
哎可惜之前没有发现论坛那么多资源。就瞎买了..太不淡定了

作者
Author: sobereva 时间: 2015-2-16 16:30
结构就直接在原有结构上自行更改就行，不用同源模建
至于拓扑文件信息，如果涉及的基团在力场里本来就有，就自行拼一下。如果比如是一个奇特的基团，找不到现有的，就截取出来，用prodrg、ATB之类生成相应的参数和拓扑信息。高斯并没法直接算。

作者
Author: kunkun 时间: 2015-2-16 18:40

sobereva 发表于 2015-2-16 16:30
结构就直接在原有结构上自行更改就行，不用同源模建
至于拓扑文件信息，如果涉及的基团在力场里本来就有， ...

谢谢sob老师，我去尝试尝试！

作者
Author: kunkun 时间: 2015-2-19 19:46

sobereva 发表于 2015-2-16 16:30
结构就直接在原有结构上自行更改就行，不用同源模建
至于拓扑文件信息，如果涉及的基团在力场里本来就有， ...

sob老师，我使用了ATB重新去计算特殊氨基酸的残基的的拓扑文件。利用itp和原来蛋白pdb的原子命名方式，修改了itp文件后，对应gromacs54a7中rtp重新编写。后来使用pdb2gmx也能生成对应的top。说明我的修改没有问题。但是我注意到了，使用ATB计算的拓扑的电荷分布（焦谷氨酸）和原先标准氨基酸（脯氨酸）的电荷分布不同，gromos貌似对五元环的电荷分配做过调整（基本C N 都为0，除了羧基C和O）。我想请问一下sob老师，对于这种情况。我应该如何重新计算优化新残基的电荷呢？ps：在ATB计算中，有几个二面角和键角的能量提示过高...这计算的结果还可信么？
还是说初步的计算结果，我需要用高斯09在进行优化和计算电荷？
但是这样计算出来的电荷是单独氨基酸的，但是我感觉rpt文件中的氨基酸电荷分布是按照形成酰胺键的二肽为基础的..

附上氨基酸结构：有一羰基不同。 (, 下载次数 Times of downloads: 85) (, 下载次数 Times of downloads: 81)

祝sob老师们羊年快乐！

作者
Author: sobereva 时间: 2015-2-19 21:47
gromos力场库里的原子电荷不是直接计算出来的，都是手工调过，一方面是为了满足charge group定义的要求，一方面也是为了使电荷有可移植性。
你可以用AM1-BCC算一下这两种情况的电荷分布，然后再对照gromos力场里脯氨酸已有的电荷来估量羰基和相邻碳的电荷应该怎么设。或者也可以直接参考gromos力场里已有的带羰基的分子，看看羰基区域电荷是什么样的，能否直接移植到你的体系上。

你先看一下对这个氨基酸与脯氨酸一致的部分，ATB给出的键、角的参数和gromos力场里的脯氨酸是否一致。然后，对于涉及到羰基的参数，把ATB分配的力常数放到gromos力场文件里搜一下，看看对应的键角/二面角项涉及的原子类型是否和当前实际情况一致，如果一致，并且实际动力学跑起来它的构象没什么异常就行了。

作者
Author: kunkun 时间: 2015-2-19 22:13
本帖最后由 kunkun 于 2015-2-19 22:14 编辑

sobereva 发表于 2015-2-19 21:47
gromos力场库里的原子电荷不是直接计算出来的，都是手工调过，一方面是为了满足charge group定义的要求，一 ...

sob老师我不是很懂该如何搜索力场内的力常数..我之前在修改rtp时，发现标准氨基酸的参数都是以ga_xxx代表.
并无具体的数值..tpr文件也无法打开.

。。手册中也没有明说，这些参数的来源是哪..

——————
发现是在bond.itp中有相应的代码啦！谢谢sob老师

作者
Author: kunkun 时间: 2015-3-2 15:13

sobereva 发表于 2015-2-19 21:47
gromos力场库里的原子电荷不是直接计算出来的，都是手工调过，一方面是为了满足charge group定义的要求，一 ...

sob老师，我在这段时间内做了许多的动力学（关于烷烃微粒对酶的激活研究），我尝试了很多参数也尝试了更加精确地控温和控压，还有模型的建立。但是稳定性太不理想了。
比如，我同一个模型经过模拟退火升温和限制性模拟控压比较稳定后。不同时间的两次重复模拟（没有任何参数的重新确定，直接重新grompp）但结果出来却相差径庭。

1 我想请问下老师，有无做多次重复的必要？（我看很多酶方面的模拟，参数设置大同小异，也和我的差不多，但也只有一条轨迹的参数，无重复）。还是说我直接跑出比较理想的轨迹就接着往下面研究，做对接？
个人感觉分子模拟的时间仅仅10几ns的，也不太可能每次都模拟出想要的蛋白变化。（特别是更加精确地口袋变化的研究）

2 面对这种稳定性不高的情况下，有什么方法可以提高我的重复性？（关乎参数等等方面的设定？）
我看一些文献平衡100ps就完事了，连压力是否稳定都不提及....（按照我的经验 100ps完全是不可能会压力收敛的啊，大分子体系）

3 我在md模拟的过程中，发现烷烃分子对酶内部的形状影响实在太大，有时候烷烃钻入口袋时结果还不错，但是烷烃不进去的时候...结果完全不能看，这种情况下，我能否现初步获得一定构象变化的蛋白，然后以这个模型人为把分子对接进去后，继续模拟完全活性的状态？这样是否太人为了？

作者
Author: sobereva 时间: 2015-3-2 16:19
分子动力学模拟就是可重现性很差。如果模拟的过程是那种凭直觉100%能出现的还好，比如一个水簇一加电场会拉长、一个配体做SMD会被拉出去等等，怎么模拟都差不多；而对于那种自己凭直觉说不好会模拟出什么结果的情况，运气成分太大了。

动力学模拟其实放在论文中的那部分模拟本身耗时不多，但是来回模拟摸索参数所花的时间可能要比这多一个数量级以上。我当年最初做MD时模拟一个酶，其实就是跑10ns而已，但翻来覆去摸索、跑了一年有余才得到想要的现象。很多文章看起来模拟从技术角度没什么难度，但自己想效仿着做立刻跑出想要的漂亮结果，并不是简单事。

1 就先顺着比较理想的那个轨迹接着做吧。或许比如跑10次能有6次想要的，但没准儿运气差，跑几次都碰上结果不好的，若来来回回跑时间都耽误了。既然文献里也没反复重复模拟，你不去重复也有理由。

2 这个很难说。可重复性关键还是取决于体系特征。真想完全可重复那就得用串行计算，也不随机生成速度之类，但这样的话每次结果精确一致，也没重复跑的意义了。
100ps是太短了，除非真是跑不动（比如早年间的模拟），否则起码也得500ps。有些体系，比如混合膜，甚至得冲着100ns跑才能完全平衡。

3 这属于说法问题。往好的地方解释，可以说是增强采样、加速化学过程、突出重点什么的，但得能清楚证明烷烃肯定能进去。

作者
Author: kunkun 时间: 2015-3-2 20:20

sobereva 发表于 2015-3-2 16:19
分子动力学模拟就是可重现性很差。如果模拟的过程是那种凭直觉100%能出现的还好，比如一个水簇一加电场会拉 ...

每次和老师交流完都觉得更有动力了！

老师我还有小疑问：
我用insert-molecules生成一个充满烷烃的box（5 5 5），然后进行升温模拟。但最后10ns后的密度都会比正常的高。进行空压后盒子本身也在发生变化。然后我以xx.gro在vmd中选出一个所需要的半径的球再按照老师生物膜的教程，把蛋白放进去。这种做法是否合理呢？因为按照很多文章所说是在md过程中让烷烃自身集聚成球，我觉得这个摸索条件不好做，同时也太随机了，经常会模拟成烷烃柱子。（我现在有点顾虑是否是我这种方法所平衡的烷烃，其实不太真实的，因为密度都比实验值偏大，特别是C8烷，密度还要高出100g/ml）

作者
Author: sobereva 时间: 2015-3-3 04:23
一开始让烷烃分布区域距离盒子边界有足够的距离，使之不和镜像碰上，就不会成为柱子。

模拟烷烃的过程本身没什么问题，和实验有偏差和往往是力场的问题。

那么放蛋白没问题，但更好的做法是先把蛋白放进去再去除与蛋白过近接触的分子，否则可能空隙大。

作者
Author: kunkun 时间: 2015-3-3 15:54

sobereva 发表于 2015-3-3 04:23
一开始让烷烃分布区域距离盒子边界有足够的距离，使之不和镜像碰上，就不会成为柱子。

模拟烷烃的过程本 ...

老师说的是gromacs力场参数问题还是我的小分子力场参数有问题？（我是用的是gromos54a7）。小分子使用的是ATB生成的。没有用ambertools和高斯拟合resp...

作者
Author: sobereva 时间: 2015-3-3 16:01

kunkun 发表于 2015-3-3 15:54
老师说的是gromacs力场参数问题还是我的小分子力场参数有问题？（我是用的是gromos54a7）。小分子使用的 ...

主要是非键作用参数，包括LJ参数和原子电荷。
原子电荷是自己或ATB产生的。LJ参数直接来自gromos力场文件。

作者
Author: kunkun 时间: 2015-3-5 17:46

sobereva 发表于 2015-3-3 16:01
主要是非键作用参数，包括LJ参数和原子电荷。
原子电荷是自己或ATB产生的。LJ参数直接来自gromos力场文 ...

谢谢sob老师，我先重复几次小体系的烷烃看看结果怎么样。
老师我有一个小问题，我从pdb数据库下载下来的晶体是具有ca离子的活性中心结构的。对于这种情况，我该如何处理pdb文件呢？我看论坛上有些人说要用高斯算参数和电荷..

作者
Author: sobereva 时间: 2015-3-5 18:44

kunkun 发表于 2015-3-5 17:46
谢谢sob老师，我先重复几次小体系的烷烃看看结果怎么样。
老师我有一个小问题，我从pdb数据库下载下来 ...

没有配位键的情况，就当金属是个普通离子处理就行了。如果其余部分没什么特殊性，该用什么参数还是什么参数。
如果发现模拟过程中离子跑了，可以再加上限制势。

作者
Author: kunkun 时间: 2015-3-6 13:47
本帖最后由 kunkun 于 2015-3-6 14:08 编辑

sobereva 发表于 2015-3-5 18:44
没有配位键的情况，就当金属是个普通离子处理就行了。如果其余部分没什么特殊性，该用什么参数还是什么参 ...

谢谢sob老师~
但是我直接pdb2gmx会提示力场参数中无ca和na的参数，对于这种附带金属离子的是否需要先把金属离子拆分出来，生成gro后再重新把坐标复制过去？（但发现坐标对应不上了）
老师所说的限制势是否指的是限制动力学中的定义xyz限制力的itp？

如果蛋白晶体内部含有水分子，那是去除好呢还是不去除好呢？我看有些文章是保留的，有些文章又clean了，众说纷纭啊。是不是保留结晶水时，我觉得md会更加稳定一些，重复性也会好点呢？（毕竟内部有水分子可能可以稳定酶的口袋结构）..虽然我一直是把水处理掉了

作者
Author: sobereva 时间: 2015-3-6 14:08
限制势是指的原子间的距离限制

水的问题无所谓，反正加水后还会覆盖那些区域。就算保留了水，之后跑起来也和后加的水无法分辨了。除非某些水的位置和作用特别特殊，那么应该留着。

作者
Author: kunkun 时间: 2015-3-7 02:26

sobereva 发表于 2015-3-6 14:08
限制势是指的原子间的距离限制

水的问题无所谓，反正加水后还会覆盖那些区域。就算保留了水，之后跑起来 ...

谢谢sob老师，刚看了看我的一个模型，2个md重复的初步结果，我发现活性位点内的水对我的构象变化有是些帮助，看来是不能去掉的。

老师说的限制原子间的距离，我看了力场里水分子的的settle限制。想尝试模仿写一个。但是这些限制都是基于在一个分子内部（水分子，甲醇分子等），若ca等金属离子并没有与蛋白形成配位键的话，那该如何编写呢？

作者
Author: sobereva 时间: 2015-3-7 07:37

kunkun 发表于 2015-3-7 02:26
谢谢sob老师，刚看了看我的一个模型，2个md重复的初步结果，我发现活性位点内的水对我的构象变化有是些帮 ...

先尝试模拟，如果Ca没自己跑掉，就不用限制势。
如果需要加限制势，把Ca和蛋白作为整个一个分子考虑。gmx里分子间没法加限制势。

作者
Author: lonemen 时间: 2018-1-11 15:14
很好的经验和tips，谢谢分享！

作者
Author: zzffzz33 时间: 2023-4-14 20:01

sobereva 发表于 2015-1-19 17:00
1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

卢老师，想请教一下您回复的第一条，我在13 13 13 的盒子里跑含80个左右原子的有机小分子，PBC盒子里填充了300个小分子，在进行动力学计算之前先进行了平衡计算，发现PBC盒子的边界处的分子散架了，分成小片段，当时没发现这个问题，在平衡计算进行之后直接进行了分子动力学计算，结果显示盒子内的分子结构正常，而边界处的分子散架，请问这个是什么原因造成的？是否影响体系模拟正常的进行？

作者
Author: sobereva 时间: 2023-4-15 00:36

zzffzz33 发表于 2023-4-14 20:01
卢老师，想请教一下您回复的第一条，我在13 13 13 的盒子里跑含80个左右原子的有机小分子，PBC盒子里填充 ...

弄清楚是真的散架还是只是分子被盒子截断了而显得散架

倘若是真正的散架，也不会只恰好出现在盒子边界

欢迎光临计算化学公社 (http://bbs.keinsci.com/)