gromacs能量最小化出现问题！！

tjuchan

刚接触gromacs，想模拟一下小分子有机物(混合物)。先拿苯练手！
遇上一些问题，想请教一下各位！！
先讲一下我使用过程(避免是操作上的失误！！)
1、得到(PRODRG)苯的pdb文件(benzene.pdb)，建立一个空盒子（empty.gro）；
2、向空盒子中添加200个苯分子  genbox -cp empty.gro -ci benzene.pdb -nmol 200 -o benzene.gro ;
3、编写top文件和mdp文件，使用的是oplsaa力场，如附件！不知道写对没！！

***这里遇到一个问题是苯它是平面分子，我的力场给了improper 二面角，但是还是观察苯上的原子不共面***

4、grompp -f  benzene.mdp -c benzene.gro -p system.top -o benzene.tpr
5、mdrun -v -deffnm benzene
6、 按理第2步得到的初始结构并不是合理的，需要先进行能量最小化步。但是得到的结果发现。能量最小过程后(超过迭代步数，能量并未收敛值规定值)，苯分子的构型完全变了，如图1,2  (前后比较！！)

yakun777

水和苯的compressibility不一样，也许能改善计算得到的rho

sobereva

xpyp 发表于 2017-7-27 22:06
求教，在作酶－配体复合物的能量最小化，在Amber中经常采用的策略是：
1) 固定蛋白质和配体，仅优化溶剂 ...

amber教程里搞得过于复杂了，虽然那种做法很保险，但一般没必要那么啰嗦
gmx里用freeze可以冻结指定部分

xpyp

sobereva 发表于 2015-8-24 02:42
没有你的itp文件，但一看2.png就知道，绝对是苯的itp文件不对，氢原子都跑到苯环里面去了。
力场设定正常 ...

求教，在作酶－配体复合物的能量最小化，在Amber中经常采用的策略是：
1) 固定蛋白质和配体，仅优化溶剂和离子；
2)固定酶，优化配体，溶剂和离子
3)固定酶的主链，优化其它；
4)放开约束，作整体优化。
而在gmx中经常只做最后一步，是什么原因？
另外 前三步带约束的优化，在gmx中如何实现呢？

tjuchan

wbn 发表于 2016-5-9 22:49
可能npt的时间不够吧，密度来不及平衡，你才做了50 ps，太短了，至少得做几 ns 吧。先npt 再 nvt 是没问 ...

谢谢，我多试几次

wbn

tjuchan 发表于 2016-5-9 10:52
非感谢您和sob老师，按照您的建议，仔细的寻找了下原因。发现是自己的大意，采用prodrg产生的不是oplsaa ...

可能npt的时间不够吧，密度来不及平衡，你才做了50 ps，太短了，至少得做几 ns 吧。先npt 再 nvt 是没问题的，但我觉得盒子里全是苯的话全程npt也没问题，不需要nvt嘛。为什么设rlist > rvdw我也不懂，我一般设它们相等的。

tjuchan

wbn 发表于 2016-4-25 03:43
opls-aa 力场里CA-HA的力常数是 307105.6， 你用了400000, 大了不少，改正后可能解决这个warning, 如果不能 ...

非感谢您和sob老师，按照您的建议，仔细的寻找了下原因。发现是自己的大意，采用prodrg产生的不是oplsaa而是gromos的力场参数。于是又下载了mktop制作苯分子的top文件，键的力场数都对应上了。
如ben.top文件所示。但是（298k + 1atm的密度）模拟值（1.004）和文献值（0.878），相差甚远。不知道是mdp参数设置错误还是，模拟步骤错误。故再次请教各位。
我的计算步骤如下[参考gromacs的官网教程]
EM(能量最小化)---NVT(100ps)----NPT(100ps)---NPTequilibrium(5ns)
，我将这几部的mdp文件贴出来大家看看。

--------------------------
--------------------------
也看文献，看看别人的模拟细节是怎样设计的，但是有一篇文献（见附件）写的稍微详细点，但是发现跟manual上介绍的有些出入。请教一下大家，这篇文献写的模拟细节有问题吗
1、“Each system was simulated for 4 ns in an NVT ensemble with an integration time step of 0.001 ps. This was preceded by a 100 ps of NPT simulation to get the correct density.”
它的模拟步骤是NPT（100ps, 达到正确密度）+NVT(4ns，平衡步)。这是不是反了
2、“The cut-off distance of 1.2 nm was used for Lennard-Jones potential. The Coulomb potential was calculated using Particle Mesh Ewald39 with a cut-off of 1.3 nm and Fourier grid spacing of 0.12. The neighbor list was updated every 0.01 ps within 1.3 nm.”，
按照他的设置，rlist>rvdw。manual上写的是VDW采用cut-off时候，rvdw≥rlist。这。。

wbn

opls-aa 力场里CA-HA的力常数是 307105.6， 你用了400000, 大了不少，改正后可能解决这个warning, 如果不能的话就用sob神的方法。还有就是二面角，在opls-aa力场里improper dihedrals 是不能取代 proper dihedrals的，得先把proper dihedral写完，然后再在成环原子上加improper dihedral。 我不知道你拓扑文件里那些参数都哪来的，反正不是opls-aa的常数。top文件的内容一定要一条一条和文献参数核对，看看有没有错漏，不能软件生成的往里一粘贴的就不管了。建议先去看一看别人的文件是怎么写的 http://www.gromacs.org/Downloads ... ighlight=topologies

sobereva

tjuchan 发表于 2016-4-24 17:49
之前学习gromacs 5.0模拟苯的问题一直没有得到解决，就一直搁在那里。最近重新捡起来。但是在NPT步骤遇到了 ...

图1中断裂问题，你在VMD里显示周围的镜像盒子，如果此时边界分子看起来不再断裂就是PBC原因，不用管它。
那个warning是因为你的步长是2fs，但是有没有对氢涉及的键进行约束，其振动频率太高，用2fs的步长容易出问题。把constraints设hbonds就行了。

tjuchan

之前学习gromacs 5.0模拟苯的问题一直没有得到解决，就一直搁在那里。最近重新捡起来。但是在NPT步骤遇到了一些问题。现在总结下，希望大神指点一下。（最初级）主要是看了gromacs官网的第一个教程和常见error中关于‘Residue 'XXX' not found in residue topology database’的问题解答。
模拟对象是298K下液态苯（有机小分子，opls-aa）的性质（密度等）
1、先是在GV下建立苯的pdb文件（附件benzene.pdb），直接用gmx pdb2gmx 好像不行，于是进行了修改（附件benzene1.pdb），同时也增加了些东西在oplsaa.ff文件下的*.rpt文件（附件aminoacids.rtp）；
2、利用packmol进行建模[这一步不是很懂，为啥我采用gmx solvate -cs benzene1.pdb -box 25.0 25.0 25.0就不行]，packmol过程的文件是input.in和benzene.pdb，当然计算过100个苯0.847g/ml下的盒子大小25A左右。这样建模完后的文件benzene2.pdb
3、gmx -f benzene2.pdb -o benzene2.gro（选择15opls-aa，7不含水力场），修改benzene.gro文件的盒子大小（最后一行的三个数都改为其中最大值）。这一步得到的topol.top文件不靠谱主要是二面角类型是3，但是苯平面应该采用improper二面角才是，于是用PRODRG产生了合适的二面角，见文件PRODRG.top。它产生的二面角类型就是2，然后将topol.top中的二面角按照PRODRG.top进行修改为benzene.itp保存到oplsaa.ff下面。
4、新建一个system.top：包含oplsaa力场和benzene.itp，同时设置分子为100个
5、网上随便抄了个minim.mdp；gmx grompp -f minim.mdp -c benzene2.gro -p system.top -o em.tpr
6、gmx mdrun -v -deffnm em； 这一步能量最小化，问题是能量才-85，但是力满足要求。看例子说能量应该很低的。
7、VMD看了一下最小化之后的em.gro构型，发现有些分子（盒子边缘）出现断裂，这是不是因为周期边界的原因呢？见图1
8、之后再进行 gmx grompp -f nvt.mdp -c em.gro -p system.top -o nvt.tpr却出现
“WARNING 1 [file system.top, line 10]:
  The bond in molecule-type BEN between atoms 1 CA1 and 7 HA1 has an
  estimated oscillational period of 9.6e-03 ps, which is less than 5 times
  the time step of 2.0e-03 ps.
  Maybe you forgot to change the constraints mdp option.”
的错误，不知道怎么搞得。希望大家指点一下。

sobereva

tjuchan 发表于 2015-8-24 20:59
非常感谢sob老师。收获颇多。
按照group的意思是假如我需要计算苯中C-C的rdf，就连续选择相同的group(只 ...

C-C的rdf是这样。
C-H的你并不需要动itp，你在index文件里定义新的group，一个group里面是所有C的编号，另一个是所有H的编号，然后计算rdf时分别选这两个group就行了。

tjuchan

sobereva 发表于 2015-8-24 19:46
你之前怎么得到的2.png？你的trr文件最后一帧和2.png明显不一样，现在的优化的trr文件是正常的。
你跑完 ...

非常感谢sob老师。收获颇多。
按照group的意思是假如我需要计算苯中C-C的rdf，就连续选择相同的group(只包含C原子)就行吧，若是计算C-H的rdf我需要提前在itp文件中将C和H分开编到不同的group吧，是这个意思吗？

sobereva

tjuchan 发表于 2015-8-24 16:31
您好，我是采用mktop_2.2.1工具直接生成的oplsaa力场的itp文件，应该不会错吧！
另外请教一下，采用g_rd ...

你之前怎么得到的2.png？你的trr文件最后一帧和2.png明显不一样，现在的优化的trr文件是正常的。
你跑完MD了么？当前在真空中高度分散的状态计算rdf没有意义。

g_rdf的参考组就是指使用哪个group里的原子作为计算rdf的参考位置（横坐标为0的位置），之后选的group就是把哪个group里的原子作为rdf的考察对象。

tjuchan

sobereva 发表于 2015-8-24 02:42
没有你的itp文件，但一看2.png就知道，绝对是苯的itp文件不对，氢原子都跑到苯环里面去了。
力场设定正常 ...

您好，我是采用mktop_2.2.1工具直接生成的oplsaa力场的itp文件，应该不会错吧！
另外请教一下，采用g_rdf做径向分布函数<g_rdf -f benzene.trr -s benzene.tpr -o rdf.xvg>需要经历选取group，如下
Select a reference group and 1 group
Group     0 (         System) has  2400 elements
Group     1 (          Other) has  2400 elements
Group     2 (            ben) has  2400 elements




其中的reference group和1 group指的是?
我没有在定义什么group。做出来的rdf数据(峰值)，好几百都有。

sobereva

没有你的itp文件，但一看2.png就知道，绝对是苯的itp文件不对，氢原子都跑到苯环里面去了。
力场设定正常、初始结构合理的话肯定不会有Too many LINCS warnings这种提示

下面是我做最小化时常用的输入文件模板，一般都在上面根据实际情况稍作修改。
define =
integrator = steep
emtol = 100.0
emstep = 0.01
nsteps = 10000
; output frequency
nstxout = 10
nstvout = 0
nstfout = 0
nstlog = 50
nstenergy = 50
nstxtcout = 0
xtc_grps = system
;
nstlist = 10
pbc = xyz
rlist = 1.2
cutoff-scheme = group
coulombtype = PME
rcoulomb = 1.2
vdwtype = cut-off
rvdw = 1.2
DispCorr = EnerPres
;
constraints = none
constraint_algorithm = LINCS
implicit_solvent = no

其实对于你的模拟体系，分子之间离得那么远，不会有不合理接触，不优化都没关系。

编译时默认就是编译出能单机并行的，你一看CPU占用率就知道。计算时默认会用机子上所有核心，人为指定要用的核数的话用-nt选项。