Gaussian可识别的GAFF和GAFF2力场文件

snljty

Gaussian软件的势能面扫描等功能非常强大，可以多种方式指定限制性优化等功能。但是起内置的分子力场种类较少，目前只有Amber，Dreiding和UFF力场，缺乏常见的描述有机小分子的GAFF力场。虽然用Gromacs（AmberTools应该也有很多类似功能，这里不提）的rerun即可完成力场下的单点能计算功能，但是对于柔性扫描（限制型优化），就显得比较麻烦，通常做法是对要被限制的变量加一个非常大的限制势，然后做收敛阈值很小的能量极小化。不过虽然原理上能量极小化和量子化学常说的几何优化很相似，终究目的不一样，而且Gromacs这方面收敛阈值也不如Gaussian条件更多。如果可以在Gaussian中使用GAFF力场进行限制型优化，也许会方便一些。考虑到GAFF和Amber力场形式完全相同，Gaussian有内置的Amber力场，原理上应该可行。为此，从AmberTools20里面找到了gaff.dat和gaff2.dat的力场文件，对照Gaussian16自带的amber.prm，写了个Python脚本将其转化为Gaussian能接受的形式。
使用的时候请自行指定原子类型，连接方式，原子电荷。将相关的力场文件在合适的位置用@引用即可。如果不了解怎么引用，可以参考卢老师博文http://sobereva.com/60 第5小节。几何优化类任务一定要有Opt=NoMicro关键词。力场的关键词要有Amber=SoftOnly。

欢迎反馈。


files.zip (106.8 KB, 下载次数 Times of downloads: 271)

2021-6-7 13:26 上传 Uploaded

点击下载
Click to download

exity

高能！！！

ginlpein

请问大佬有在Gaussian上用GAFF力场做柔扫的输入文件范例吗？我采用您上传的gaff2.prm做柔扫，但是一直报L101错误，卡在L101这里无法顺利运行。
我输入文件格式：

1. %chk=XXXXX.chk
   
2. #p geom=connectivity opt=(modredundant,NoMicro) amber=SoftOnly
   
3. 【空一行】
   
4. XXXXXX
   
5. 【空一行】
   
6. 0 1
   
7. 【原子坐标】
   
8. 【空一行】
   
9. 【原子键接信息】
   
10. 【空一行】
    
11. D 72 73 76 93 S  30 3.00000
    
12. 【空一行】
    
13. @./gaff2.prm

复制代码

输出的log尾部：

1. Generating MM parameters.
   
2. Read MM parameter file:
   
3. Define C3         1
   
4. End of file in RdPar.
   
5. Error termination via Lnk1e in /share/apps/g16/A03/l101.exe at Thu Nov  3 22:21:16 2022.
   
6. Job cpu time:       0 days  0 hours  0 minutes 13.2 seconds.
   
7. Elapsed time:       0 days  0 hours  0 minutes  3.4 seconds.
   
8. File lengths (MBytes):  RWF=      6 Int=      0 D2E=      0 Chk=      1 Scr=      1
   

复制代码

输出的log文件开头：

1. Entering Gaussian System, Link 0=g16
   
2. Input=RTPEcage_opt.com
   
3. Output=RTPEcage_opt.log
   
4. AtFile(1): ./gaff2.prm
   
5.             ! Generated by FF_Amber2Gaussian.py
   
6.             !
   
7.             ! Non-bonded interaction
   
8.             !
   
9.             NonBon 3 1 -10 0 0.000 0.000 0.500 0.000 0.000 -1.200
   
10.             !
    
11.             ! Stretches
    
12.             !
    
13.             HrmStr1 ow hw  553.0  0.9572
    
14.             HrmStr1 hw hw  553.0  1.5136
    
15.             HrmStr1 br br   81.0  2.5420
    
16.             HrmStr1 br c1  227.0  1.7870
    

复制代码

程序刚运行进入L101时的输出：

1. (Enter /share/apps/g16/A03/l101.exe)
   
2. -------------------
   
3. XXXXXX
   
4. -------------------
   
5. Symbolic Z-matrix:
   
6. Charge =  0 Multiplicity = 1
   
7. 【原子笛卡尔坐标】
   
8. 【空一行】
   
9. The following ModRedundant input section has been read:
   
10. D 72 73 76 93 S  30 3.00000
    
11. 【空一行】
    
12. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
13. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
14. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
15. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
16. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
17. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
18. ITRead=  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    
19. 【省略】

复制代码

从输出文件看似乎力场是正常读取并加载了，原子笛卡尔坐标也有正确读入，但是程序始终卡在L101出错，不知道是哪里格式出问题了。
有考虑可能是原子键接信息这块导致读取出问题，但是力场计算键接信息不可少，也不可能删掉来避错。
请问大佬是否有成功的输入案例借我参考一下？