用Psi4计算RESP电荷
1. RESP电荷计算
有很多不同方案能把分子的电荷分布投射到原子上，例如Mulliken、Löwdin以及NPA电荷等。有一类基于拟合静电势的电荷，如CHELPG、Merz-Kollman (MK)和RESP电荷。RESP电荷因在AMBER和GAFF力场中的使用而闻名。
能直接计算RESP电荷的量化计算软件并不多。往往需要额外的波函数分析工具去辅助计算RESP电荷。例如可以通过Gaussian/GAMESS与R.E.D.(The RESP and ESP charge Derive)的联合使用计算，或用Multiwfn分析波函数信息计算RESP，或者用Amber提供的工具进行计算。
Psi4量子化学计算包提供了RESP计算模块，让我们能从头到尾用一个软件完成整个的RESP电荷计算。这个计算模块提供用户定义约束等价原子电荷的功能以及对多个构象进行加权处理。作者采用Connolly分子表面生成格点进行计算。可以进行两步优化，即先不开启对称性约束以及小的约束强度进行对极性部分进行拟合，然后第二步拟合加上对称性限制，拟合CH3以及CH 2的电荷分布。
如果用户有Psi4使用经验或者会Python，很容易就能进行计算。没有相关经验或者只会Python的简单使用，可以用我的计算模板，稍微修改就能快速计算。
2. Psi4的使用（PsiAPI为例）
Psi4的安装可以参考公众号之前的文章《PSI4程序安装及运行》，建议下载预编译版本或者用Docker镜像运行程序。我用的版本是Psi4conda-1.4rc3-py39-Linux-x86_64.sh。另外，Psi4还支持wsl2版以及有windows版，建议大家试试。
Psi4有两种使用模式，一种是类似大多数常规程序的可执行文件运行输入文件控制方法（Psithon），另一种是写个Python脚本与Psi4交互的方式（PsiAPI）。我目前只会第二种，而且第二种方式更加灵活。
2.1 参数设置
参数分为两大部分，一个是全局参数设置，一个是RESP计算的参数设置。
全局参数设置控制整体量化计算是开壳层还是闭壳层。通过psi4.set_options方法指定。例如
psi4.set_options({'reference': 'rks', # 闭壳层Kohn-Sham计算,
'scf_type' : 'direct'}) # 不采用density-fitting，直接计算电子积分更多设置请参照手册。RESP参数设置如下：
options =
{'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0], #不同壳层格点生成设置(建议默认)
'VDW_POINT_DENSITY' : 1.0, #格点密度，1点每Å^2(建议默认)
'RESP_A' : 0.0005, #第一步参数，越小约束越弱(建议默认)
'RESP_B' : 0.1, #约束参数(建议默认)
'RESTRAINT' : True, #RESP计算必须开启(建议默认)
'IHFREE' : False, #是否约束氢，False为对氢进行约束。
'CONSTRAINT_GROUP' : [[4,2],[8,6],[11,3],[12,7]], #原子[4,2],[8,6]... 定义为等价原子
'WEIGHT' : [1, 1, 1], #三个构象，构象的权重
'METHOD_ESP' : "b3lyp", #计算级别
'BASIS_ESP' : "cc-pvtz" #基组
}用户需根据体系指定开壳层或闭壳层计算（'scf_type'），以及指定等价原子（'CONSTRAINT_GROUP'），其余的建议默认。
2.2 分子结构设置
以有三个构象的CF 3CH2 OH为例，三个构象通过如下方式指定：
# Initialize three different conformations of ethanol
geometry = """
0 1
C 0.00000000 0.00000000 0.00000000
C 1.48805540 -0.00728176 0.39653260
O 2.04971655 1.37648153 0.25604810
H 3.06429978 1.37151670 0.52641124
H 1.58679428 -0.33618761 1.43102358
H 2.03441010 -0.68906454 -0.25521028
F -0.40814044 -1.00553466 0.10208540
F -0.54635470 0.68178278 0.65174288
F -0.09873888 0.32890585 -1.03449097
"""
mol1 = psi4.geometry(geometry)
mol1.update_geometry()
mol1.set_name('conformer1')
geometry = """
0 1
C 0.00000000 0.00000000 0.00000000
C 1.48013500 -0.00724300 0.39442200
O 2.00696300 1.29224100 0.26232800
H 2.91547900 1.25572900 0.50972300
H 1.61500700 -0.32678000 1.45587700
H 2.07197500 -0.68695100 -0.26493400
F -0.32500012 1.02293415 -0.30034094
F -0.18892141 -0.68463906 -0.85893815
F -0.64257065 -0.32709111 0.84987482
"""
mol2 = psi4.geometry(geometry)
mol2.update_geometry()
mol2.set_name('conformer2')
geometry = """
0 1
C 0.00000000 0.00000000 0.00000000
C 1.48805540 -0.00728176 0.39653260
O 2.04971655 1.37648153 0.25604810
H 3.06429978 1.37151670 0.52641124
H 1.58679428 -0.33618761 1.43102358
H 2.03441010 -0.68906454 -0.25521028
F -0.27127997 -0.97245518 -0.41089913
F -0.60950912 0.20545158 0.87999334
F -0.17244493 0.77215758 -0.74975690
"""
mol3 = psi4.geometry(geometry)
mol3.update_geometry()
mol3.set_name('conformer3')
molecules = [mol1, mol2, mol3]分子自旋多重度、电荷与坐标指定方式与Gaussian的输入文件约定相同。用户可以根据上述模板自行修改或添加、减少构象。
2.3 运行与结果得出
计算的整体脚本如下，脚本修改自Psi4的测试文件。这个脚本写了完整的RESP两步计算：
import psi4
import resp
# Initialize three different conformations of ethanol
geometry = """
0 1
C -0.2753000000 -2.2228000000 0.0000000000
C 1.0407000000 -1.7594000000 0.0000000000
C -1.0775000000 0.0564000000 0.0000000000
C -1.3343000000 -1.3149000000 0.0000000000
H -0.4776000000 -3.3036000000 0.0000000000
H 1.8753000000 -2.4754000000 0.0000000000
H -1.9127000000 0.7719000000 0.0000000000
H -2.3715000000 -1.6800000000 0.0000000000
C 0.5221000000 2.0333000000 0.0000000000
H 1.5322000000 2.3862000000 0.0000000000
C 1.2975000000 -0.3886000000 0.0000000000
H 2.3254000000 -0.0916000000 0.0000000000
C 0.2381000000 0.5197000000 0.0000000000
O -0.4298000000 2.8564000000 0.0000000000
"""
psi4.set_options({'reference': 'rhf',
'scf_type' : 'direct'})
mol1 = psi4.geometry(geometry)
mol1.update_geometry()
mol1.set_name('conformer1')
molecules = [mol1]
# Specify options
options = {'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0],
'VDW_POINT_DENSITY' : 1.0,
'RESP_A' : 0.0005,
'RESP_B' : 0.1,
'RESTRAINT' : True,
'IHFREE' : False,
'CONSTRAINT_GROUP' : [[4,2],[8,6],[11,3],[12,7]]
# Call for first stage fit
charges1 = resp.resp(molecules, options)
print("Restrained Electrostatic Potential Charges")
print(charges1[1])
options['RESP_A'] = 0.001
resp.set_stage2_constraint(molecules[0], charges1[1], options)
options['grid'] = []
options['esp'] = []
# Add constraint for atoms fixed in second stage fit
for mol in range(len(molecules)):
options['grid'].append('%i_%s_grid.dat' %(mol+1, molecules[mol].name()))
options['esp'].append('%i_%s_grid_esp.dat' %(mol+1, molecules[mol].name()))
# Call for second stage fit
charges2 = resp.resp(molecules, options)
print("\nStage Two\n")
print("RESP Charges")
print(charges2[1])执行 python test.py > out 计算结果会被重定向到out文件中，在out最后一行，可以看到RESP电荷：
Stage Two
RESP Charges
[ 0.5935773 0.06749319 -0.6303203 0.40224024 0.0254598 0.0254598
-0.16130334 -0.16130334 -0.16130334]程序开发者测试显示，Psi4的RESP计算结果与R.E.D. 计算结果误差小于10−4 。
3. 更多计算模板
用户的计算是非常简单的，只需要将自己的体系套入到模板中，运行即可快速得到结果。脚本可以在以下地址中下载，
https://github.com/chenxin199261/QM-Templates/tree/main/Psi4-RESP
我为大家准备了开壳层、闭壳层计算以及多个构象计算的模板。
使用Psi4计算RESP不要忘记引用RESP模块以及Psi4本体程序。
4. 参考资料
https://github.com/cdsgroup/resp
https://onlinelibrary.wiley.com/doi/10.1002/qua.26035