计算化学公社

标题: 使用sobtop建立聚合物的itp以后，如何分配RESP电荷到每个原子 [打印本页]

作者
Author: 牧生 时间: 2022-2-23 20:13
标题: 使用sobtop建立聚合物的itp以后，如何分配RESP电荷到每个原子
本帖最后由牧生于 2022-2-23 20:24 编辑

sobtop得到的聚合物itp，是没有分配原子的电荷的。以一个嵌段聚合物，聚合度均为50的聚乙二醇-聚甲基丙烯酸甲酯（图省事，简写为PEG50-PMMA50）为例，请教如何去添加正确的RESP电荷。聚合物结构式如图所示：
(, 下载次数 Times of downloads: 59)
该聚合物结构明确，且元素很常见，所以直接全部使用了GAFF，几秒钟就得到itp文件。

第一步：用M$画一个聚合度为3的PEG，以羟基封端，然后用orca进行Optimization，溶剂为水，关键词为
! B97-3c opt noautostart miniprint nopop
%maxcore  1000
%pal nprocs 4 end
%cpcm
smd true
SMDsolvent "water"
end
* xyz 0 1
这里有一个疑问，Opt以后，可以得到一个gbw文件，用这个也可以计算RESP电荷，但是潜意识觉得哪里有点不对，觉得应该是在溶剂中优化结构，然后再在溶剂情况下算单点能才是正确步骤。

第二步：然后用Multiwfn打开得到的.out文件，输出单点能计算的inp文件，关键词为如下
! B97-3c noautostart miniprint nopop
%maxcore  1000
%pal nprocs 4 end
%cpcm
smd true
SMDsolvent "water"
end
* xyz 0 1

第三步: 用Multiwfn打开得到的gbw文件，计算RESP电荷

RESP电荷结果如下：

O -6.664604 0.825225 0.189018  -0.6980158077
C -5.477920 1.588142 0.371967 0.2473893021
C -4.295225 0.715090 0.015145 0.3569874671
H -5.397759 1.920393 1.413500 0.0047525467
H -5.494722 2.478285 -0.265392 0.0047525467
H -4.348984 0.421843 -1.031528 0.0135556298
H -4.305610 -0.193427 0.622525 0.0135556298
O -3.063094 1.420507 0.175531  -0.7845293930
H -2.946468 1.589966 1.118905 0.4615227089
O -10.196692 1.401893 0.590626  -0.7271621912
C -9.022682 0.658007 0.293611 0.3170818095
C -7.836117 1.567119 0.500735 0.4869132643
H -9.050058 0.300767 -0.740861  -0.0117338890
H -8.948130 -0.214271 0.950746  -0.0117338890
H -7.802343 1.914250 1.538416  -0.0307126182
H -7.914520 2.445765 -0.147586  -0.0307126182          //红色字体这一段是中间的-O-CH2-CH2-，即为PEG中的大量重复单元，该片段的总电荷-0.00806
O -13.732324 0.686135 0.681205  -0.8291114464
C -12.558512 1.499695 0.766593 0.2349083906
C -11.374066 0.615490 0.464297 0.4784937587
H -12.454913 1.927369 1.767439 0.0142139403
H -12.602644 2.320946 0.047084 0.0142139403
H -14.491552 1.266042 0.811138 0.5094995843
H -11.454346 0.209750 -0.549446  -0.0170643331
H -11.344733 -0.226980 1.163332  -0.0170643331

计划用中间那个-O-CH2-CH2-的电荷值作为PEG50中重复单元的电荷值，手动调一下，把-0.00806分配到O和C上，使得该片段的总电荷为0。那么，PEG50中的大量重复单元就可以使用这组数据。

第四步：以相同的方法，画一个PMMA3的分子式，用相同的方法获取RESP电荷，
PMMA3的RESP电荷结果如下：

C -8.861314 0.275272 1.163147  -0.0693266705
C -7.332534 0.311744 1.147430  -0.1846717283
H -9.218124 0.843117 0.306200 0.0764443119
H -9.193342 0.811810 2.047335 0.0764443119
C -6.768534 0.190134 -0.244847 0.9636955882
H -6.922187 -0.498651 1.742891 0.0708831338
O -7.191673 0.777122 -1.228215  -0.7050741735
O -5.681967 -0.600338 -0.278629  -0.4264116191
C -5.018055 -0.726665 -1.554696  -0.0548554677
H -5.690297 -1.162177 -2.288968 0.1008191691
H -4.177101 -1.386948 -1.379985 0.1008191691
H -4.669237 0.242607 -1.899867 0.1008191691
C -6.827244 1.636393 1.726102  -0.2957661952
H -7.247169 2.481581 1.182931 0.0917242579
H -5.741591 1.705158 1.681465 0.0917242579
H -7.130560 1.716078 2.767916 0.0917242579
C -11.048981 -0.967867 0.840688  -0.2678499309
C -9.533049 -1.100328 1.128247  -0.0903302656
H -11.480126 -1.966063 0.900847 0.1402297498
H -11.140978 -0.658899 -0.198235 0.1402297498
C -9.345348 -1.903402 2.410286 0.8258256270
O -9.779067 -3.035695 2.561236  -0.6688996571
O -8.645221 -1.271259 3.358174  -0.2064877541
C -8.462863 -1.975847 4.603762  -0.2709968095
H -7.915364 -2.899927 4.441293 0.1446951737
H -9.425639 -2.190556 5.059461 0.1446951737
H -7.889866 -1.305516 5.233003 0.1446951737
C -8.959718 -1.969526 -0.001608  -0.5063934677
H -9.016076 -1.429472 -0.943726 0.1483215751
H -9.528410 -2.889963 -0.103196 0.1483215751
H -7.922622 -2.240604 0.182445 0.1483215751
H -13.693554 -1.301914 1.489498 0.1346688467          //这一段是中间的甲基丙烯酸甲酯，总电荷为0.109046
C -13.395731 -0.285305 1.239397  -0.4461449506
C -11.944797 -0.052532 1.699728 0.0187367893
H -13.502195 -0.151313 0.165772 0.1346688467
H -14.077234 0.400199 1.741032 0.1346688467
C -11.664799 1.425361 1.471289 0.8602932354
O -11.524056 2.256836 2.354194  -0.6843512889
O -11.648373 1.748151 0.169211  -0.3181579126
C -11.395344 3.131251 -0.149543  -0.2840809859
H -10.426251 3.437593 0.235162 0.1567653284
H -12.174489 3.766500 0.262717 0.1567653284
H -11.401867 3.185181 -1.231600 0.1567653284
C -11.849694 -0.364355 3.186649  -0.3470808292
H -10.885594 -0.091866 3.604471 0.1077047189
H -12.010237 -1.427743 3.348119 0.1077047189
H -12.612596 0.175639 3.743026 0.1077047189
计划用中间那段作为PMMA50中大量重复单元的原子电荷，且把0.109046平均分配到其余各个原子，使得该类片段总电荷为0。  如果某个片段带电为+1或者-1，那么也用相似的方法处理，使得该片段的总电荷为+1或-1。。

第五步：PEG与PMMA的连接处，再自己手动建立一个PEG2-PMMA2去计算电荷，再手动复制到PEG50-PMMA50相应的itp中连接处的位置。端基处也用相同方法处理，并且我觉得，端基不严格正确也没关系，毕竟绝大多数重复单员都是正确的。这样，就可以得到一个电荷正确的PEG50-PMMA50的itp文件。

请教一下，我的操作有没有不合理的地方，或者是否有多余或遗漏，该怎么改正呢？

作者
Author: sobereva 时间: 2022-2-24 11:20
要求不高的情况，可以直接让Sobtop指认MMFF94电荷，最为省事

对于你描述的过程也没有明显问题
但不建议自己手动调电荷使得单元净电荷成为整数，应当在Multiwfn算RESP电荷的时候施加电荷约束条件直接令单元的净电荷为整数，这样对静电势重现性比你手动调会更好，下文有获得氨基酸残基RESP电荷的例子可参考
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）

另外，对于这种没有局部显离子性的体系，优化时没有绝对必要带溶剂模型。ORCA里带溶剂模型后耗时会高不少，而且亦有可能由于数值噪音影响收敛（虽然ORCA 5.0在这方面改进了不少）。

作者
Author: 牧生 时间: 2022-2-25 17:52

sobereva 发表于 2022-2-24 11:20
要求不高的情况，可以直接让Sobtop指认MMFF94电荷，最为省事

对于你描述的过程也没有明显问题

谢谢解答，已经会用约束让片段的电荷为0了。
现在还有两个问题请教一下：
第一、
如果用orca计算波函数，那么，直接就Multiwfn载入pdb文件，直接使用如下设置得到的波函数质量，再计算出来的RESP电荷，
-1 Toggle employing implicit solvation model, current: No
0 Select task, current: Single point
1 B97-3c 或者用 1b r2SCAN-3c
和调用gaussian的懒人脚本计算得到的RESP电荷，哪个质量高一些呢？？因为看懒人脚本里面，是先优化，再在溶剂下计算单点能，所以始终没想通。
限于本人就是手残党，是量化外行，还没参加过量化培训，社长能否写一个调用orca计算RESP电荷的懒人脚本呢？

第二、

对于这种没有局部显离子性的体系，优化时没有绝对必要带溶剂模型。

复制代码

那么，如果是带电荷的片段，如丙烯酸的，那么是不是要带溶剂模型，如果不带溶剂模型那又会怎样呢？

作者
Author: sobereva 时间: 2022-2-25 22:18

牧生发表于 2022-2-25 17:52
谢谢解答，已经会用约束让片段的电荷为0了。
现在还有两个问题请教一下：
第一、

1 懒人脚本的更高一些
不过r2SCAN-3c的静电势质量也已经挺不错了
ORCA的懒人RESP脚本我早就打算写，一直没时间

2 应当带。局部显离子性的体系的构型受溶剂效应影响可能较大

作者
Author: 高处裹棉被 时间: 2022-5-4 23:45
楼主好，请问对于头尾两个片段的RESP原子电荷是如何计算得到的？

比如用PEG3约束中间片段电荷为0后，计算得到结果后，将羟基封端的那一个片段原子电荷就作为整个聚合物头部片段的RESP原子电荷吗？

谢谢！

作者
Author: 牧生 时间: 2022-5-5 07:32
本帖最后由牧生于 2022-5-5 07:34 编辑

高处裹棉被发表于 2022-5-4 23:45
楼主好，请问对于头尾两个片段的RESP原子电荷是如何计算得到的？

比如用PEG3约束中间片段电荷为0后，计 ...

要约束端基也为相应的电荷，比如PEG的端基电荷是0
http://sobereva.com/soft/Sobtop/ 例5：产生聚合度为40的氯丁二烯聚合物的拓扑文件

操作方法见
RESP拟合静电势电荷的原理以及在Multiwfn中的计算
http://sobereva.com/441（http://bbs.keinsci.com/thread-10880-1-1.html）

作者
Author: 一起搞科研 时间: 2023-4-18 22:47
本帖最后由一起搞科研于 2023-4-18 22:53 编辑

牧生发表于 2022-5-5 07:32
要约束端基也为相应的电荷，比如PEG的端基电荷是0
http://sobereva.com/soft/Sobtop/ 例5：产生聚合 ...

楼主，请问你在套用中间重复单元的时候有没有遇到问题啊？我使用gaussview生成.mol2文件后用sobtop构建top文件，但是套用中间单元电子的时候非常困难，麻烦。因为中间单元的原子序号是不规律的，而我的约束的RESP电荷的文件的原子序号是规律的。请问您有遇到过这样的问题吗？

作者
Author: 牧生 时间: 2023-4-19 08:37
本帖最后由牧生于 2023-4-19 08:45 编辑

一起搞科研发表于 2023-4-18 22:47
楼主，请问你在套用中间重复单元的时候有没有遇到问题啊？我使用gaussview生成.mol2文件后用sobtop构建to ...

按照我以前用过的经验，使用gv建立了分子模型，如果中间单元的原子序号是不规律的，就不会出现“约束的RESP电荷的文件的原子序号是规律的”
sobtop以及Multiwfn不会改变gv建立分子的原子编号以及类型等。但是，如果建立的分子模型不合理，倒是有可能让sobtop或者Multiwfn辨认不准确，或者编号发生变化。举例：一个分子上有好几个羟基或者羧基或者羰基，如果建模时，某些原子靠的太近，比如羰基和临近的羧基太近，就会导致羰基氧原子和羧基氢原子“成键”。把这个结构载入sobtop时，就会有一个原子类型变成问号。

作者
Author: 一起搞科研 时间: 2023-7-31 15:46

牧生发表于 2023-4-19 08:37
按照我以前用过的经验，使用gv建立了分子模型，如果中间单元的原子序号是不规律的，就不会出现“约束的RE ...

谢谢大佬回复

作者
Author: a-Student 时间: 2023-8-4 00:32
本帖最后由 a-Student 于 2023-8-4 10:09 编辑

sobereva 发表于 2022-2-24 11:20
要求不高的情况，可以直接让Sobtop指认MMFF94电荷，最为省事

对于你描述的过程也没有明显问题

老师打扰您一下，我构建5个单体聚合，想要计算RESP电荷根据您的博文步骤，我把中间的单体施加电荷约束为0，但是计算出来中间单体的总电荷不是0，请问老师可能是什么原因呢？文件已上传，fchk文件太大了所以上传网盘了，辛苦老师指教，谢谢老师。
[已解决，谢谢老师]

作者
Author: sobereva 时间: 2023-8-4 00:34

a-Student 发表于 2023-8-4 00:32
老师打扰您一下，我构建5个单体聚合，想要计算RESP电荷根据您的博文步骤，我把中间的单体施加电荷约束为0 ...

8成是你在Multiwfn里选择了标准两步式RESP拟合，此时约束不会对两步都生效，界面上的文字和屏幕上的提示都明确体现了

作者
Author: wwwv 时间: 2025-7-16 09:06
请问我使用Multiwfn+orca优化和计算好单体的电荷了，怎么使用sobtop构建多聚体啊

作者
Author: flymetothemoon 时间: halfhour ago

a-Student 发表于 2023-8-4 00:32
老师打扰您一下，我构建5个单体聚合，想要计算RESP电荷根据您的博文步骤，我把中间的单体施加电荷约束为0 ...

您好，您这个问题是因为采取了两步拟合吗？您是将聚合物的每个重复单元的电荷都设置成了0吗？

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