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

sobereva

没搞懂你的意思。

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

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

kunkun

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

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

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

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

sobereva

1nm也无妨。
反复试，总能找到合适的设定。

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

kunkun

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

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

kunkun

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

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



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轨迹有感。

sobereva

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

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

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

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

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

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

kunkun

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

谢谢sob老师！
同感啊，心情和结果挂钩..做了30多个模型，跑出了2个..也算欣慰了
最后一天在实验室，希望做出点结果好回家过年！

kunkun

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老师，新年快乐！

sobereva

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

kunkun

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

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

sobereva

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

kunkun

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

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

kunkun

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文件中的氨基酸电荷分布是按照形成酰胺键的二肽为基础的..

附上氨基酸结构：有一羰基不同。



祝sob老师们羊年快乐！

sobereva

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

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

kunkun

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

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


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