计算化学公社

标题: 模拟出的分子结构全乱了 [打印本页]
作者
Author: sljm    时间: 2023-10-29 23:03
标题: 模拟出的分子结构全乱了
本帖最后由 sljm 于 2023-11-1 23:14 编辑

如图所示，旁边的较小分子和右下角两个分子模拟的，完了结果如图1：MD结果1，分子跑乱了， (, 下载次数 Times of downloads: 8) ，但是跑的过程前半段分子结构还是正常的如图2：MD中间状态2 (, 下载次数 Times of downloads: 8) 这是运行全程报错过的信息 (, 下载次数 Times of downloads: 8)                              warning1是电荷计算0.001，我把它忽略了，1AKI分子原子数53即分子1（结果pdb截图），unk即pdb截图中的2[attach]77814[/attach]，相关文件已上传，请问是哪里出错了呢？或者我需要跑更多的步数运行吗还是截取最后一步正常结构状态即可

作者
Author: laoman    时间: 2023-10-30 01:15
没有给log信息，也没有给grompp的warnning信息。个人经验觉得这是grompp生成MD的tpr那一步，输入的PDB/gro里的小分子顺序和itp（力场）文件里的对应不上。我这里有一个简单的sort PDB原子顺序的脚本：
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-                                                                                                                                                     
  3. import sys                                                                                                                                                                  
  4. # Usage: python script.py LIG_GMX.itp LIG_old.pdb LIG_new.pdb                                                                                                               
  5. COM = [line for line in open(sys.argv[2], 'r') if line.split()[0] in ['ATOM', 'HETATM']]                                                                                    
  6. RES = COM[0].split()[3]                                                                                                                                                     
  7. ITP = [[x for x in line.split()] for line in open(sys.argv[1], 'r') if RES in line and len(line.split()) > 8]                                                               
  8. OUT = open(sys.argv[3], 'w')                                                                                                                                                
  9.                                                                                                                                                                            
  10. ITP_ATM = [x[4] for x in ITP]                                                                                                                                               
  11. COM_ATM = [x.split()[2] for x in COM]                                                                                                                                       
  12.                                                                                                                                                                            
  13. if len(ITP_ATM) != len(COM_ATM):
  14.     print('Number of atoms not match!')
  15.     exit(1)

  17. NEW = [COM[COM_ATM.index(ITP_ATM[x])] for x in range(0, len(ITP_ATM))]
  18. OUT.write(''.join([x for x in NEW]))
复制代码


lig_old.pdb 可以这样获得：
  1. gmx editconf -f com_solv_ions.gro -o com_solv_ions.pdb
  2. grep "$LIG_residue_name" com_solv_ions.pdb > lig_old.pdb
复制代码

然后用lig_new.pdb里的内容去替换com_solv_ions.pdb里的对应内容。
最后重新跑grompp：
  1. gmx grompp -f md.mdp -c com_solv_ions.pdb -r com_solv_ions.pdb -n index.ndx -p topol.top -o md.tpr
复制代码


如果grompp那一步没有任何warnnings请忽略我的回答
作者
Author: sljm    时间: 2023-11-1 23:19
laoman 发表于 2023-10-30 01:15
没有给log信息，也没有给grompp的warnning信息。个人经验觉得这是grompp生成MD的tpr那一步，输入的PDB/gro ...

谢谢解答，我检查并重新运行了程序，报错信息我重新上传了，但是是电荷的报错信息部位0，计算为0.001，我给忽略了，相关的文件和信息我又上传了一遍，结果还是很奇怪
作者
Author: laoman    时间: 2023-11-2 00:41
本帖最后由 laoman 于 2023-11-2 00:49 编辑
sljm 发表于 2023-11-1 23:19
谢谢解答，我检查并重新运行了程序，报错信息我重新上传了，但是是电荷的报错信息部位0，计算为0.001，我 ...

嗯，你这warrning并没有提及原子顺序不对应的情况。提取PDB应该用trjconv：gmx trjconv -s name.tpr -o md.pdb -f md.gro/xtc -pbc mol -n index.ndx
你的情况估计是PBC没有去掉，用trjconv处理的时候加个-pbc mol试试。
有时候这个小体系的PBC反而不好处理。
作者
Author: sljm    时间: 2023-11-3 20:21
[/img]
作者
Author: sljm    时间: 2023-11-3 20:33
laoman 发表于 2023-11-2 00:41
嗯，你这warrning并没有提及原子顺序不对应的情况。提取PDB应该用trjconv：gmx trjconv -s name.tpr -o m ...

嗯嗯，我后来发现可能是因为我constraints在文件中设置为了h-bonds，由于分子跑的太开氢键被限制所以键被拉伸了，我把这个设置成all-bonds后分子结构没有再发生这种情况，不过两个分子还是没有形成稳定的构象，这在我的实验中似乎有一点不太正常，这是我后来的md.pdb最后一帧截图，和RMSD，所以我想请问的是有没有办法把两个分子限制在一定的范围内进行模拟而不会跑的太散，或者设置更长的步数可能会有结果，或者别的方法，目前跑的步数5000000，10ns (, 下载次数 Times of downloads: 9) (, 下载次数 Times of downloads: 8)
作者
Author: laoman    时间: 2023-11-3 23:07
sljm 发表于 2023-11-3 20:33
嗯嗯，我后来发现可能是因为我constraints在文件中设置为了h-bonds，由于分子跑的太开氢键被限制所以键被 ...

”我把这个设置成all-bonds后分子结构没有再发生这种情况，不过两个分子还是没有形成稳定的构象“。这就足够了。至于模拟的结果不符合实验结果，可能的因素有很多，比如选用的力场，小分子的电荷设置，初始速度等等，MD有时候就是个玄学的东西。你这个体系不大，我建议可以直接用Gaussian的DFT+色散校正去算相互作用。Gaussian里面怎么算弱相互作用，你可以搜一下sob老师的帖子（http://sobereva.com/）



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3