使用sobtop建立聚合物的itp以后，如何分配RESP电荷到每个原子

牧生

sobtop得到的聚合物itp，是没有分配原子的电荷的。以一个嵌段聚合物，聚合度均为50的聚乙二醇-聚甲基丙烯酸甲酯 （图省事，简写为PEG50-PMMA50）为例，请教如何去添加正确的RESP电荷。聚合物结构式如图所示：

该聚合物结构明确，且元素很常见，所以直接全部使用了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文件。


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

wwwv

请问我使用Multiwfn+orca优化和计算好单体的电荷了，怎么使用sobtop构建多聚体啊

sobereva

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

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

a-Student

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

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

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

一起搞科研

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

谢谢大佬回复

牧生

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

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

一起搞科研

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

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

牧生

高处裹棉被 发表于 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）

高处裹棉被

楼主好，请问对于头尾两个片段的RESP原子电荷是如何计算得到的？

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

谢谢！

sobereva

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

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

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

牧生

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电荷的懒人脚本呢？

第二、

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

复制代码

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

sobereva

要求不高的情况，可以直接让Sobtop指认MMFF94电荷，最为省事

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

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