有关两相体系模拟的问题请教

kunkun

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



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

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

2015-1-19 15:56 上传 Uploaded

点击下载
Click to download
mdp

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中的应用是否合适，因为我看有人发过相关的文章）

sobereva

1 不清楚你的意思，什么叫跑散。
是否是因为PBC造成的，看这些片段是否都在边界就清楚了。

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

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


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


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

kunkun

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

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

kunkun

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..

sobereva

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

kunkun

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的存在是为了缓冲区间？

sobereva

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设得大一些也就基本计算精确了。

kunkun

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

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

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

kunkun

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 一起使用.........难道是软件之间对此概念的设计不同么？


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

kunkun

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

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

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

sobereva

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

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

kunkun

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

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

kunkun

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的选项，但是他的解释是分子会分开输出.我就没敢用）

kunkun

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

找到解决办法了..

kunkun

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

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

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

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

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


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