关于使用Gaussian对PMMA低聚物进行结构优化的问题

光荣道彭于晏

卢老师您好。
1，我最近在学习如何使用pdb2gmx建立线性聚合物的拓扑文件。我想先计算出PMMA中间重复单元以及封端两个残基的RESP电荷，于是先用Chemdraw3D绘制了结构为ABBBC的五聚体的.mol2结构，然后用Gaussview保存为.gjf。
关键词写为“# opt b3lyp/6-311g** scrf=(solvent=ethanol)”，关键词是否合理？溶剂是否需要设置？如果需要，要设置成何种溶剂？（这种小分子的聚集态一般都是液态的）



2，我本来计划先用acpype.py先处理一下五聚体的结构，然后以产生的目录中的五聚体_NEW.pdb的结构为基础，再进行Gaussian的结构优化（电脑性能有限，希望通过这种方式减少一点计算时间）。但是使用VMD查看该结构，发现图中标出的位置成键了，三个原子分别是C-H-O，按说这应该是一个氢键；同时，GV中发现不光有C-H-O，还有C-H-H-C。这种情况是否是脚本运行出错？此结构是否可以作为Gaussian结构优化的基础？

光荣道彭于晏

sobereva 发表于 2022-1-8 06:31
量化计算之前，必须肉眼检查结构合理性，别上来就算，否则可能得到毫无意义的结构
如果存在不合理的近距离 ...

卢老师您好，非常感谢您提供的解答和帮助。我按照您提供的教程，一步步摸索，终于完成了RESP2电荷的计算。

一、计算过程（写这个主要是想给像我一样的纯小白用户一个参考，另外也是想让卢老师给把把关）

1. 绘制残基顺序为ABC的MMA三聚体。
之前使用ChemDraw和ChemDraw3D绘制，导致最终出现乱连键的情况，所以尝试直接使用GaussView直接绘制。但是结构优化之后还是出现了C-O-O-C这种键。根据卢老师的提示，这种情况应该和所用软件没有关系，而是要在GaussView中仔细检查分子，看看是不是有距离过近的原子，通过二面角工具手动调整分子的构象，使之有一个相对合理的初始结构，原子之间尽量疏离一些。保存为.gjf文件。
（对于一点也不会Gaussian和GaussView的同学，可以先看一下Gaussian入门教程，纯小白向）


2. 在.gjf文件中编辑关键词。
这一步非常重要，卢老师教程里面是将三个任务（结构优化、气相下的计算和液相下的计算）写到一个文件里面了，而我分成了三个任务分别跑（其实没有这个必要）。我在初期尝试计算的时候出现了报错，没有意识到需要去对应的.out文件里查找详细的报错原因（对于.out文件的简单解读，在Gaussian入门教程视频里面也有，终端提示的报错信息完全无法确定问题所在），就分成了三个任务分别跑，尝试确定错在了哪个任务。

结构优化：

1. %chk=/run/media/root/HDD2_Linux-Win/MD_HDD2/MMA-1-RESP2(0.5)/opt.chk
   
2. # B3LYP/TZVP em=GD3BJ opt

复制代码

气相下：

1. %oldchk=/run/media/root/HDD2_Linux-Win/MD_HDD2/MMA-1-RESP2(0.5)/opt.chk
   
2. %chk=/run/media/root/HDD2_Linux-Win/MD_HDD2/MMA-1-RESP2(0.5)/SP_gas.chk
   
3. # B3LYP/def2TZVP em=GD3BJ geom=allcheck

复制代码

液相下：

1. %oldchk=/run/media/root/HDD2_Linux-Win/MD_HDD2/MMA-1-RESP2(0.5)/opt.chk
   
2. %chk=/run/media/root/HDD2_Linux-Win/MD_HDD2/MMA-1-RESP2(0.5)/SP_solv.chk
   
3. # B3LYP/def2TZVP em=GD3BJ scrf(read,solvent=generic) geom=allcheck
   
4. 
   
5. eps=4.0
   
6. epsinf=1.9

复制代码

关键词的设置方法，参照楼上卢老师提供的链接。

3. 调用Gaussian进行计算

1. g16 < 文件名.gjf > 文件名.out          （我自己安装的是Gaussian16，就用g16）

复制代码

通过检查三个任务的.out文件，证明所有的计算都已完成。

结构优化：

1.                  Item                 Value           Threshold    Converged?
   
2. Maximum   Force                0.000024      0.000450       YES
   
3. RMS          Force                 0.000003      0.000300       YES
   
4. Maximum  Displacement      0.001669      0.001800       YES
   
5. RMS          Displacement      0.000353      0.001200       YES
   
6. Predicted change in Energy=-1.061396D-08
   
7. Optimization completed.
   
8.     -- Stationary point found.          （出现4个yes就证明完成了结构优化）

复制代码

1. Job cpu time:       2 days 10 hours 28 minutes 27.0 seconds.
   
2. Elapsed time:       0 days  3 hours 41 minutes 18.7 seconds.   （翻到最下面能看到任务持续的时间，第二个时间是任务持续的实际时长。i7-10700K，设置Gaussian使用16核，47个原子的结构，在B3LYP/TZVP下，任务持续不到四小时，完全可以接受）

复制代码

经过检查.chk文件，发现没有乱连键的情况，而且分子的姿态更为舒展，原子之间的距离也更加疏离，应该是得到了正确的优化结构。


气相下：

1. Job cpu time:       0 days  3 hours 21 minutes 23.0 seconds.
   
2. Elapsed time:       0 days  0 hours 12 minutes 45.6 seconds.     （在较好的优化结果基础上，气相下的计算很快就能完成）

复制代码

液相下：

1. Job cpu time:       0 days  4 hours 27 minutes 16.7 seconds.
   
2. Elapsed time:       0 days  0 hours 16 minutes 57.3 seconds.      （液相下的计算同样很快）

复制代码

4. 计算RESP2(0.5)电荷

（参照教程在centOS下正确安装Multiwfn       https://www.bilibili.com/video/av41402462/）

把SP_gas.chk和SP_solv.chk都转换为fch文件：

1. formchk SP_gas.chk

复制代码

1. formchk SP_solv.chk

复制代码

将/sob/Multiwfn_3.8_dev_bin_Linux/examples/RESP/calcRESP.sh 拷贝到当前目录   （这里偷懒使用了脚本来计算RESP2(0.5)电荷）

1. ./calcRESP.sh SP_gas.fch SP_solv.fch   （这行代码是直接复制卢老师帖子的，但是我这里报错了，报错原因是因为上一步的格式转换，把.chk转成了.fchk，所以这里也必须是.fchk，不然脚本找不到文件）

复制代码

1. ./calcRESP.sh SP_gas.fchk SP_solv.fchk  （这一行代码就能正确使用脚本并调用Multiwfn来计算了）

复制代码

“运行过后，RESP电荷会输出到当前目录下的与输入文件同名但是后缀为chg的文件中，而RESP2电荷会输出为RESP2.chg。chg是Multiwfn的私有的记录原子电荷的格式，其中2~4列是来自fch文件里的坐标，最后一列是算出来的电荷值。”

5. 确定残基的原子和原子电荷，把数据和分子结构对应起来

在CentOS中，使用GaussView打开opt.chk文件（我的win系统下，GaussView6.0搭配Gaussian09打不开该文件），右击—view—labels，让程序显示原子序号；然后观察分子结构，找到封端的残基和中间的残基，并通过分别右击原子——选择原子的方式，讲中间的残基标注出来（选中之后的原子会变成黄色）

将RESP2.chg文件中原子顺序与GaussView打开的结构进行对照，发现两者是完全对应的。（这里学过量化的同学应该可以直接能确定，我是自己手动验证的）

将RESP2.chg数据拷贝到表格文件，并进行分列，对照GaussView，把其中三个残基的原子分别标注出来。到这里就完成了RESP电荷计算。

二、几个问题需要再请教卢老师

1. RESP2电荷的计算包含了液相下的计算，对于我要计算的任务，需要自定义溶液参数。而甲基丙烯酸甲酯三聚体是液体，PMMA是固态，分子结构上两者只有重复单元数这一点区别，那么，在自定义溶液的eps和epsinf时，要参照甲基丙烯酸甲酯单体溶液的参数（三聚体的参数应该很难找），还是要参照PMMA固体的相关参数？

2. 另外，时间关系，我只在一篇中文期刊上查到了PMMA的介电常数是4.0，根据您的帖子，把epsinf设成了1.9（“大部分溶液的epsinf在这附近”），这个参数是否合理？

3. 还有就是，对于MMA这种体系，应该不需要考虑液相非极性部分贡献吧？

光荣道彭于晏

sobereva 发表于 2022-1-8 06:31
量化计算之前，必须肉眼检查结构合理性，别上来就算，否则可能得到毫无意义的结构
如果存在不合理的近距离 ...

谢谢卢老师的耐心解答，我先研究一下您发来的这几个链接的内容。

sobereva

量化计算之前，必须肉眼检查结构合理性，别上来就算，否则可能得到毫无意义的结构
如果存在不合理的近距离接触，必须手动调节二面角转开。或者基于正确的连接关系用分子力场做优化。只有当严重不合理接触都消除掉，肉眼看上去像个正经的分子结构，才能用量化方法去优化

先用acpype.py没有丝毫意义。你要搞清楚每个程序是以什么方式做什么。acpype处理过程中结构会被ambertools里的SQM在AM1级别下优化，这个级别早就过时了，而且SQM速度又慢。倘若你要快速做个预优化，应当用主流的半经验方法如PM7（G16、MOPAC都支持）或者GFN-xtb（xtb程序支持）做预优化。对于你这种体系，做个分子力场级别的预优化也足够了。

仔细看此文了解优化用的方法问题，不带色散校正显然不能用
谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的
http://sobereva.com/557（http://bbs.keinsci.com/thread-17899-1-1.html）

构造模型体系算RESP电荷用三个重复单元就够了，没必要用5个。封端的单元太多，不仅优化耗时高，还可能优化完了结构卷起来给中间那个关注的单元碍事。

当前不是在乙醇下算的，显然不适合solvent=ethanol。自定义溶剂，设成这种聚合物的实际的介电常数。怎么自定义下文说了
谈谈隐式溶剂模型下溶解自由能和体系自由能的计算
http://sobereva.com/327（http://bbs.keinsci.com/thread-3345-1-1.html）
优化带不带溶剂模型无所谓，但给Multiwfn算RESP电荷用的波函数必须是在考虑溶剂模型下产生的，因为溶剂模型对电荷分布有显著的极化作用，下文说了
RESP2原子电荷的思想以及在Multiwfn中的计算
http://sobereva.com/531（http://bbs.keinsci.com/thread-16190-1-1.html）

wzkchem5

这些多余的键说明Chem3D画出的结构不合理。画完以后应该先用分子力学之类的方法初步优化一下，更好的做法是退火一下，然后再用高斯做计算。
计算理论级别必须加色散校正，且优化完了必须做freq证明是局部极小值点结构，此外应该没什么问题