计算化学公社

标题: 计算RESP原子电荷的超级懒人脚本（基于schrodinger的Jaguar模块和Multiwfn程序） [打印本页]

作者
Author: wangzhehyd 时间: 2022-1-9 17:22
标题: 计算RESP原子电荷的超级懒人脚本（基于schrodinger的Jaguar模块和Multiwfn程序）
本帖最后由 wangzhehyd 于 2022-1-9 17:24 编辑

药物设计方向的同仁会比较熟悉Schrodinger软件包，此贴参考社长的《计算RESP原子电荷的超级懒人脚本（一行命令就算出结果）》 http://bbs.keinsci.com/thread-12858-1-1.html，利用Schrodinger中的Jaguar模块和Multiwfn程序中的RESP电荷计算模块编写了一个实现RESP原子电荷计算的懒人脚本，可用于批量计算（无需手动指定分子的电荷和多重自旋度，貌似程序可以预判）。

module load schrodinger multiwfn
先载入Schrodinger和Multiwfn环境变量，集群中安装了Lmod，直接module load相应程序即可

run ./run_jaguar_resp.py h2o.sdf 1
run为schrodinger程序中的命令，run_jaguar_resp.py为所编写的懒人脚本，h2o.sdf为输入的分子的结构，1为计算时调用的cpu核心数

纯属量化小白瞎折腾，烦请懂量化的专家提改进和优化的意见。

(, 下载次数 Times of downloads: 38)

作者
Author: 牧生 时间: 2022-1-10 09:41
本帖最后由牧生于 2022-1-10 09:59 编辑

帖子里面有一句话，(无需手动指定分子的电荷和多重自旋度，貌似程序可以预判）。
那么，预判的准确性高吗？我有时候需要用orca计算一个小体系的单点能，比如几个带电有机小分子的聚集体，总电荷倒是容易看出来，但作为懒人，我都是直接用0,1,2,3去试自旋多重度，只要报错，那就肯定不是；哪个值能算成功，就认为是它。

我的方法肯定是不合理的，那么，能用楼主帖子中的方法，傻瓜式的获得自旋多重度吗？

作者
Author: chengb917 时间: 2022-1-21 10:41
好东西，想测试一下看看。

作者
Author: sobereva 时间: 2022-1-22 16:15

牧生发表于 2022-1-10 09:41
帖子里面有一句话，(无需手动指定分子的电荷和多重自旋度，貌似程序可以预判）。
那么，预判的准确性高吗 ...

只有成键方式普通的有机体系才可能傻瓜式地判断，因为自旋多重度几乎总是越低越稳定
但碰上过渡金属体系、双/多自由基体系、卡宾等电子结构复杂的情况，那就必须自己判断了

作者
Author: wzkchem5 时间: 2022-1-22 18:01

牧生发表于 2022-1-10 02:41
帖子里面有一句话，(无需手动指定分子的电荷和多重自旋度，貌似程序可以预判）。
那么，预判的准确性高吗 ...

这样只能试出自旋多重度是奇数还是偶数。而且试成功的自旋多重度可能不止一个，你遇到成功的就停了，万一后面的自旋多重度还有能算成功的，那你仍然可能判断错。
正确的懒人方法是：在所有能算成功的自旋多重度里，取算出来能量最低的。而且要注意以下事项：
（1）只有刚进程序马上报错说自旋多重度和电子数不匹配，才属于算失败。如果因为SCF不收敛等其他原因报错，需要尝试各种方法让SCF收敛以后才能判断，而不能简单认为这个自旋多重度一定是错的。参见解决SCF不收敛的方法http://sobereva.com/61
（2）如果体系可能有反铁磁耦合，那么对于每个自旋多重度都需要准备合适的初猜（对于简单的情况可以用guess=mix，复杂的情况需要具体例子具体分析），否则直接无脑设自旋多重度算出来的能量可能是错的，比如实际应该单重态比三重态能量低，但是如果没有加guess=mix，可能会得到单重态比三重态能量高的错误结论。如果不知道怎么判断一个体系是否有反铁磁耦合，必须自己先去查阅反铁磁耦合的定义，找几个现成的有反铁磁耦合的体系熟悉一下，再回来判断自旋多重度。
（3）自旋多重度必定大于等于1，所以应该从1开始试，而不是从0开始。当然从0开始也没有错，因为程序也会告诉你0是错的，只不过你试0的时间浪费了而已。
（4）并不是所有分子都可以试到3就停的，具体试到多少可以停仍然要从化学常识来判断，靠画共振式判断这个体系最多有多少个单电子，然后才能决定试到几可以停。我见过有人算含两个高自旋Fe(III)的体系，以为试1、3、5、7就够了，结果正确答案是11。
（5）并不是所有的纯有机体系都是自旋多重度越低越稳定，甚至类药分子里也有反例。比如假如有人想用带两个spin label的药物分子做MD模拟，需要算RESP电荷，此时盲目按自旋多重度越低越稳定就有可能是错的，正确的自旋多重度很可能是3。

作者
Author: chengb917 时间: 2022-1-24 10:07
下载不了，可以发一份，测试一下看看吗？

作者
Author: aixin 时间: 2025-1-5 10:16
用楼主脚本测试了一个有机小分子，resp 和 resp2 均有原子电荷>1 和 <-1
cat 1777_resp.chg
F -16.009054 8.537424  -29.171269  -0.6933852855
C -15.623020 7.326559  -29.601532 1.2059632384
N -15.485375 7.197841  -30.898400  -0.5062942099
C -15.425672 6.323640  -28.653633  -0.6273608486
C -15.083170 6.005971  -31.376005  -0.1327848649
C -15.018125 5.074998  -29.138002 0.2215844336
C -14.815317 4.936106  -30.521005  -0.3433023738
C -15.009838 5.824346  -32.882739 0.8893311796
C -14.875629 3.893586  -28.252249  -0.1956090744
C -16.270730 5.099911  -33.377279 0.4608323405
C -14.826603 7.136561  -33.644224 1.9156643127
N -14.805002 2.725433  -28.860844 0.1097200873
C -14.857615 3.979922  -26.818158 0.6716523672
O -16.261208 3.772862  -32.852907  -0.3792132026
N -14.759229 1.665215  -28.071315  -0.2565124187
N -14.791816 2.915699  -26.047090  -1.0292461308
C -17.394781 3.064073  -32.743551 0.3802956802
C -14.756430 1.725012  -26.654949 1.4508558368

1777_resp2.chg
F    -16.009054    8.537424 -29.171269 -0.6750822545
C    -15.623020    7.326559 -29.601532 1.1698875981
N    -15.485375    7.197841 -30.898400 -0.4527855218
C    -15.425672    6.323640 -28.653633 -0.6138207991
C    -15.083170    6.005971 -31.376005 -0.1615281552
C    -15.018125    5.074998 -29.138002 0.2077543522
C    -14.815317    4.936106 -30.521005 -0.3273757903
C    -15.009838    5.824346 -32.882739 0.9053909883
C    -14.875629    3.893586 -28.252249 -0.1966494256
C    -16.270730    5.099911 -33.377279 0.4531702066
C    -14.826603    7.136561 -33.644224 1.9148668037
N    -14.805002    2.725433 -28.860844 0.1137295771
C    -14.857615    3.979922 -26.818158 0.6349897477

作者
Author: SiqiLee 时间: 2025-1-5 11:30

aixin 发表于 2025-1-5 10:16
用楼主脚本测试了一个有机小分子，resp 和 resp2 均有原子电荷>1 和

你有没有发现，你的这个小分子，它没有加H啊

作者
Author: aixin 时间: 2025-1-5 15:40
本帖最后由 aixin 于 2025-1-5 15:42 编辑

SiqiLee 发表于 2025-1-5 11:30
你有没有发现，你的这个小分子，它没有加H啊

因为原子电荷大于1和小于-1，主要是针对这个问题，所以没有贴出整个分子的结果。氢原子肯定是有的，下面是部分氢原子和重原子的电荷结果：
cat 1777_resp2.chg
F    -16.009054    8.537424 -29.171269 -0.6750822545
C    -15.623020    7.326559 -29.601532 1.1698875981
N    -15.485375    7.197841 -30.898400 -0.4527855218
C    -15.425672    6.323640 -28.653633 -0.6138207991
C    -15.083170    6.005971 -31.376005 -0.1615281552
C    -15.018125    5.074998 -29.138002 0.2077543522
C    -14.815317    4.936106 -30.521005 -0.3273757903
C    -15.009838    5.824346 -32.882739 0.9053909883
。。。
C    -17.306746    1.857405 -31.973157 -0.0911510951
C    -18.617121    3.413284 -33.293700 -0.4430155412
N    -14.723694    0.502528 -26.148064 -1.0951632485
C    -14.714964 -0.336787 -27.225606 0.5627178889
C    -14.773195 -0.103340 -29.867329 1.0871030849
C    -18.508397    1.110580 -31.770725 0.4460911422
C    -16.096726    1.417067 -31.383899 -0.4246152398
C    -19.729949    2.571547 -33.049838 0.6472354212
。。。
C    -18.443518 -0.047944 -30.948024 0.1630665748
C    -17.258697 -0.431341 -30.354065 -0.3984421119
H    -15.616753    6.530293 -27.602035 0.3201281415
H    -14.483551    3.985168 -30.925978 0.1979588924
H    -14.158932    5.157939 -33.097861 -0.4818599624
H    -16.284762    5.059773 -34.478848 -0.1511807718
H    -17.171277    5.635500 -33.035675 -0.1511807718
H    -14.700796    6.940270 -34.719139 -0.6178592201
H    -15.694702    7.798443 -33.510347 -0.6178592201
。。。
H    -15.188253    1.989100 -31.563971 0.2466600086
H    -20.693817    2.850755 -33.492582 -0.1809060502
H    -19.364635 -0.615367 -30.794853 -0.0714150315
H    -17.235254 -1.313744 -29.708447 0.1734336985

作者
Author: SiqiLee 时间: 2025-1-5 16:21

aixin 发表于 2025-1-5 15:40
因为原子电荷大于1和小于-1，主要是针对这个问题，所以没有贴出整个分子的结果。氢原子肯定是有的，下面 ...

嗷嗷，那不知道这个结果与gaussian或者ocra的计算结果差异大吗？如果比较之后差异比较大，也不排除是开壳层导致自旋多重度判断错了？

欢迎光临计算化学公社 (http://bbs.keinsci.com/)