关于为CP2K的MM模拟生成PDB文件和PSF文件的方法

gingaorb

各位小哥哥，各位大佬们，请教一个问题呀~
是这样的，我最近在模拟的体系中，包含64个水分子和一个丙酮酸分子，然后我先用CP2K的QS模块，跑AIMD，得到平衡构型的.xyz格式文件，接下来我想用MM方法继续跑动力学，于是我用VMD将该XYZ文件转换为PDB文件，参考了站长大大在CP2K培训讲义里MM部分的实例4：水中的氨基酸动力学模拟。其中的力场文件是在CHARMM22原始文件par_all22_prot.prm的基础上加入了对水TIP3P的参数，其中包括，OT-HT成键项，HT-OT-HT键角项、OT和HT的LJ参数也就是非成键项。然后我根据par_all22_prot.prm力场中的原子类型手动更改了丙酮酸分子的原子类型，并利用LigParGen在线网站，使用1.14*CM1A-LBCC电荷模型，为丙酮酸分子生成了电荷，然后让CP2K自动GENERATE CREATE_MOLECULES，这样就不需要PSF文件，但是却出现了报错，
FORCEFIELD| Missing Bond              (   C,  CC)                              
FORCEFIELD| Missing Bond              (  OC,   H)                              
FORCEFIELD| Missing Bond              (  OT,  OT)                              
FORCEFIELD| Missing Bond              (  HT,  HT)                              
FORCEFIELD| Missing Angle             ( CT3,   C,  CC)                        
FORCEFIELD| Missing Angle             (   C,  CC,  OC)                        
FORCEFIELD| Missing Angle             (  CC,   C,   O)                        
FORCEFIELD| Missing Angle             (  CC,  OC,   H)                        
FORCEFIELD| Missing Angle             (  OT,  OT,  HT)                        
FORCEFIELD| Missing Angle             (  OT,  HT,  HT)                        
FORCEFIELD| Missing Angle             (  OT,  OT,  OT)

首先Missing Bond              (   C,  CC) ，代表缺失了丙酮酸分子的羰基C和羧基CC之间的键参数，是不是意味着对于丙酮酸，我不能使用par_all22_prot.prm力场，而应该改用par_all36_cgenff.prm力场并相应的修改原子类型？

第二个问题是，对于一些小分子化合物，有没有自动的方法，基于特定的力场，例如par_all22_prot.prm力场或par_all36_cgenff.prm力场生成相应的PDB文件？还是说只能自己根据力场，手动修改每个原子的原子类型呢？


第三个问题是，对于一些小分子化合物，有没有比较好的方法去生成它的原子电荷呢？这里我用的是LigParGen在线网站，此外我看到，http://bbs.keinsci.com/forum.php ... 5&highlight=psf中，站长大大提到，机子里装了Gaussian和免费的Multiwfn的话，用下文的傻瓜式脚本直接就能通过一行命令基于pdb文件得到动力学模拟常用的RESP电荷、非常好的RESP2电荷以及很适合OPLS-AA的1.2*CM5电荷
计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）
http://sobereva.com/476（http://bbs.keinsci.com/thread-12858-1-1.html）
RESP2原子电荷的思想以及在Multiwfn中的计算
http://sobereva.com/531（http://bbs.keinsci.com/thread-16190-1-1.html）
计算适用于OPLS-AA力场做模拟的1.2*CM5原子电荷的懒人脚本
http://sobereva.com/585（http://bbs.keinsci.com/thread-21462-1-1.html）

没有Gaussian也无妨，用免费的ORCA量子化学程序也可以优化分子结构并且产生molden格式的波函数文件，可以作为Multiwfn的输入文件算原子电荷

什么电荷最合适，取决于力场。到这里我突然有个疑问，是否可以用Gaussian计算的NBO或Mulliken电荷呢？如果不行是为什么呢？

第四个问题是，对于其它报错，例如 Missing Bond              (  OT,  OT)、Missing Bond              (  HT,  HT)、Missing Angle             (  OT,  OT,  HT) 等等这些，本来就不存在呀，我感觉是不是它产生分子时，判断键接关系错了？对于这种情况，站长大大提到“本来没有的bond/angle如果出现了，说明自动猜的拓扑信息不对，应当给psf文件确保识别的连接关系正确”。那么有没有什么好的方法，可以通过结构文件来产生PSF文件呢？根据我查阅的资料，我看到VMD中的autopsf模块可以做这件事，但不幸的事，我在使用该方法时，遇到以下报错，并且无法解决，我不知道该报错的原因是什么？（我将PDB文件附在下面，如果有大佬会，希望能教教我这个小笨蛋(ૢ˃ꌂ˂ૢ)）。

然后呢，站长大大也提供了一种方法，可以用sobtop（http://sobereva.com/soft/Sobtop）产生gromacs的拓扑文件，然后parmED转成psf，分两步走。想请问一下，各位大佬，还有没有其它方法可以解决以上问题呢？希望大家可以不吝赐教哇！万分感谢大家~

sobereva

只有NPA电荷，没有NBO电荷，仔细看
计算化学中的一些常见不良写法和用词
http://sobereva.com/298（http://bbs.keinsci.com/thread-1358-1-1.html）

NPA和Mulliken电荷对静电势重现性太烂，不可能用于分子力场。仔细看《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）以及《RESP拟合静电势电荷的原理以及在Multiwfn中的计算》（http://sobereva.com/441）里的讨论。


VMD里输入以下命令可以产生只有拓扑信息而没有具体参数和原子电荷的psf文件
topo guessbonds
topo guessangles
topo guessdihedrals
topo guessimpropers
animate write psf example.psf