sljm 发表于 2023-11-3 20:33
嗯嗯，我后来发现可能是因为我constraints在文件中设置为了h-bonds，由于分子跑的太开氢键被限制所以键被 ...

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

[/img]

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反而不好处理。

laoman 发表于 2023-10-30 01:15
没有给log信息，也没有给grompp的warnning信息。个人经验觉得这是grompp生成MD的tpr那一步，输入的PDB/gro ...

谢谢解答，我检查并重新运行了程序，报错信息我重新上传了，但是是电荷的报错信息部位0，计算为0.001，我给忽略了，相关的文件和信息我又上传了一遍，结果还是很奇怪

没有给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)
    
16. 
    
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请忽略我的回答