拟合RESP电荷、力场参数

1.Gaussian计算

需要准备的是小分子gjf文件（mol，sdf等，GV打开也可保存为gjf）

1.1传统gjf文件

下面是老版的，可以用，也可以用新版的，往后翻。
gjf文件加开头结尾，最终类似下面：
  • %mem=4GB
    • 

  • %nproc=4
    • 

  • # opt hf/6-31g(d) geom=connectivity pop=mk iop(6/33=2,6/42=6,6/50=1) 
    • 


    • 

  • 3rx4_ligand
    • 


    • 

  • 0 1
    • 

  •  C(PDBName=C4,ResName=SFI,ResNum=317_A)    -8.34200000    9.07000000   17.71000000
    • 

  •  C(PDBName=C5,ResName=SFI,ResNum=317_A)    -8.95900000    8.14900000   18.68000000
    • 

  •  C(PDBName=C6,ResName=SFI,ResNum=317_A)    -5.62400000    6.49800000   18.05000000 
    • 

  •  
    • 

  •  1 15 1.0
    • 

  •  2 23 1.0 24 1.0
    • 

  •  3 5 1.0 7 2.0 10 1.0
    • 


    • 


    • 

  •  antechamber-ini.esp
    • 


    • 

  •  antechamber.esp
    • 

  •  
    • 

  •  
0 1：电荷、自旋多重度，按实际情况手动修改，自旋多重度一般都为1。
opt：几何优化。非必需，但是一般都加。个人为加的原因是，RESP电荷充分考虑了分子柔性，因此在一定空间内与构象变化关系不大，优化可以将其他手段得到的结构（晶体/动力学等）进一步优化到构象空间内与之最相近的能量极小值点。
hf/6-31g(d)：计算方法/级别。若需要溶剂，添加scrf=(solvent=xx)。考虑该参数时务必阅读4.。
geom：确保原子之间键连正确。加不加都行，理论上不影响结果，不加下面的数字键连都要删掉。
iop(6/33=2)：进行RESP Fitting并输出到Gaussian的.log文件。
IOp(6/42)：Density of points per unit area in esp fit. 0：Default (1) N：Points per unit area。
IOp(6/50)：Whether to write Antechamber file during ESP charge fitting.
pop=mk：Produce charges fit to the electrostatic potential at points selected according to the Merz-Singh-Kollman scheme 

1.2更新的gjf文件

  • %mem=4GB
    • 

  • %nproc=4
    • 

  • # opt b3lyp/6-311g(d,p) pop=mk iop(6/33=2,6/42=6) em=gd3bj
    • 


    • 

  • 3rx4_ligand
    • 


    • 

  • 0 1
    • 

  •  C(PDBName=C4,ResName=SFI,ResNum=317_A)    -8.34200000    9.07000000   17.71000000
    • 

  •  C(PDBName=C5,ResName=SFI,ResNum=317_A)    -8.95900000    8.14900000   18.68000000
    • 

  •  C(PDBName=C6,ResName=SFI,ResNum=317_A)    -5.62400000    6.49800000   18.05000000 
    • 

  •  
    • 

  •  
其实原子后面的括号里的也可以去掉。不影响结果
em=gd3bj：色散校正。对于大部分没有分子内弱相互作用（氢键，卤键等）不加影响不大，但总的来说推荐加上。
注意：这里使用的是b3lyp/6-311g(d,p)计算级别是有争议的，传统的做法是使用HF/6-31G(d)真空。具体参考4.。

1.3提交Gaussian计算

  • nohup g16 3RX4_ligand.gjf &
如果是g09就写g09，最好用g16

2.手动提取参数

2.1提取小分子的电荷

  • antechamber -i *.log -fi gout -o lig.prep -fo prepi -c resp 
–nc -1
nc后面跟小分子的formal charge，但是貌似不用这个参数也可以。
antechamber生成的以ANTECHAMBER开头的文件中的原子顺序与gaussian输入中的gjf中的原子顺序一致，而生成的NEWPDB.PDB文件和.prep文件中的原子顺序与gjf文件中原子顺序不一致，因此将ANTECHAMBER_PREP.AC的电荷信息和原子类型输出即可。
  • awk '$1=="ATOM" {print $3, toupper($10) "-" $9}' ANTECHAMBER_PREP.AC
得到：
原子名称 原子类型-电荷
N1 N3--0.297108
O1 O--0.703333
C1 C3--0.031802
复制粘贴保存为resp.txt文件，文件最后最好不要有空行，后续会用到。

2.2.提取力场信息

  • parmchk  -i lig.prep -f prepi -o forcefiled.frcmod -a y
生成的ligand.frcmod文件中有力场参数信息（可用于分子力学/动力学模拟）。
分别输入下面的命令：
  • awk '/NONBON/, /NR/' *mod | awk '{print "VDW", $1, $2, $3}'
    • 

  • awk '/^BOND/, /ANGLE/' *mod | awk '{print "HrmStr1", $1, $2, $3, $4}' | tr '-' ' '
    • 

  • awk '/^ANGLE/, /DIHE/' *mod | awk '{print "HrmBnd1", $1, $2, $3, $4, $5}' | tr '-' ' '
复制粘贴保存以上命令产生的分子力场信息，其中可能会产生一些只有VDW或者其他没有数字的行，不需要。后续会用到。

3.自动提取参数

参考《量子化学-ONIOM自动化脚本》

4.计算级别和溶剂的考虑

传统的是在真空中使用HF/6-31G(d)，是因为此级别会高估偶极矩10%-15%，正好等效反映出水的极化作用。此外，蛋白的分子力场参数拟合时，也是在这个级别下进行的，所以保持一致比较好（之前的文章基本都是这么做的）。但是，但这种高估效果并不稳定。近几年的文章使用b3lyp/6-311g(d,p)在隐式溶剂模型中计算的在不断增多。
卢天认为，传统级别在部分体系会造成较大的误差，可以使用更高级的方法进行精确计算（结合隐式溶剂模型）但是一旦结合溶剂，计算值将随着溶剂极性增大而越来越偏离实验值。比如，后续模拟是在水溶液中进行，这里也结合水模型计算，将导致RESP电荷异常，最终结果异常（因为水的极性太大了）。
综上，个人认为使用b3lyp/6-311g(d,p)。
a.如果用于ONIOM，建议使用真空。因为ONIOM就是在真空中进行的，而且对于大部分体系而言，量子化学优化几何结构时，真空中的构象和溶液中的相差不大。
b.如果用于已经考虑溶剂极化的力场，不需要使用溶剂，因为力场已经考虑这方面的信息，不需要通过原子电荷体现溶剂化效应。

5.参考