动手尝试18碳环的优化与振动分析出现的问题

不想飞的猫头鹰

    最近看卢天老师关于18碳环的博文感觉很有趣，从http://sobereva.com/515《谈谈18碳环的几何结构和电子结构》、http://sobereva.com/584《我对一篇存在大量严重错误的J.Mol.Model.期刊上的18碳环研究文章的comment》、http://sobereva.com/674《8字形双环分子对18碳环的独特吸附行为的量子化学、波函数分析与分子动力学研究》以及社长写的关于泛函与基组的博文http://sobereva.com/272《简谈量子化学计算中DFT泛函的选择》、http://sobereva.com/336《谈谈量子化学中基组的选择》中了解到，对于18碳环这样的大共轭体系，选用wb97xd这样的泛函是比较稳健的，而且社长推荐用6-311g*、def-TZVP及以上级别的基组，得到的无虚频的可靠结构应该是D9h对称性、键长交替、键角均匀160度，而6-31g*这样的不够大的基组会有问题，但是Pople基组构建得不够好，我好奇能否用def2-SVP计算出可靠的结构，另外想试试在这种已经可以当成标准结构的初入坐标输入下，6-31g*会产生什么结果，就把社长在http://sobereva.com/515提供的已经用wb97xd结合def-TZVP计算好的fchk文件下载后转换成了新的gjf文件，关键词#P wb97xd def2svp opt freq以及#P wb97xd 6-31g(d) opt freq。另外， 由于gaussian16a03的wb97xd泛函的色散校正是DFT-D2版本，我直觉上以为用orca 5.0.4中的DFT-D3版本的wb97x-d3应该会同样稳健，而且注意到可以结合 就用orca也算一下试试，关键词! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop，并使用卢老师的OfakeG工具转换orca输出文件方便gview观看有无虚频、优化出的分子结构。使用g16得到的log文件，以及使用orca的fake.out文件也上传了，感谢各位网友老师指点迷津。

    计算过程和结果如下：
wB97X-D3 TZVP得到的结构没有虚频，而且在gview里发现确实是长短键交替，但是键角却不是160度，总是有的大一丁点有的小一丁点，tools-point group里默认标准判断出的对称性是Cs。
wB97XD 6-31g*得到的结构有5个虚频，也是长短键交替，键角也不是精确的160度但是十分接近，对称性是D9h。
wB97XD def2-SVP得到的结构有1个虚频，也是长短键交替，键角也不是精确的160度但是十分接近，对称性是D9h。
然后换成6-311g*算了下，来检测出来确实是基组的问题，关键词#P wb97xd 6-311g(d) opt freq
wB97XD 6-311g*得到的结构没有虚频，也是长短键交替，键角也不是精确的160度但是十分接近，对称性是D9h。
根据卢老师的博文http://sobereva.com/278《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》，怀疑可能是在当前基组级别下优化不到理想的D9h对称性，而当前使用的是卢老师提供的满足高对称性的结构，我就按虚频模式方向略微调了下结构重新计算，终于得到了没有虚频的结构：
wB97XD 6-31g*得到的结构没有虚频，也是长短键交替，键角也不是160度，对称性也是Cs。
wB97XD def2-SVP得到的结构没有虚频，也是长短键交替，键角也不是160度，对称性也是Cs。
这说明基组确实不够大，即使换成构建得不错的def2-SVP也不能得到D9h对称性。但是很奇怪，并没有像http://sobereva.com/584中所说的那样对称性变为C9h。

按理说orca的泛函和基组都已经不错了，结果却没有优化出预期的结构，我就把基组调大到def2-TZVP,关键词! wB97X-D3 def2-TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop,但是
wB97X-D3 def2-TZVP得到的结构也没有虚频，也是长短键交替，键角也不是160度但是十分接近，对称性也是Cs。
最终尝试了wB97X-D3 def2-QZVP得到的结构没有虚频，长短键交替，键角也不是160度但是十分接近，对称性也是Cs。
我查看orca 5.0.4软件手册（加上目录的总共页数的第73页）发现wB97X-V使用DFT-D4色散校正的泛函支持解析hessian，这个泛函在一众测评中表现十分不错，就也试了下，关键词! wB97X-D4 TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop。由于wB97X-V使用VV10色散校正时不支持解析hessian就没尝试。
wB97X-D4 TZVP得到的结构没有虚频，也是长短键交替，键角也不是160度，对称性是C1。
    问题如下：
1.orca为什么没有产生正确的对称性的结构？是不是与orca在优化时没有对称性的启用有关？orca的计算级别也足够高而且带的色散版本更好，但是没有一个产生了正确的结构，是不是说，如果已经知道一个结构本该具有高对称性，使用orca这种没有启用对称性的程序进行计算是不稳健的？
2.由于怀疑是不是程序启用对称性算法的问题，g16中我用wB97XD 6-311g*加了nosymm重新计算了一下，仍然没有虚频，也是长短键交替，键角也不是160度，对称性也是Cs，这时候我就想不通，对称性算法的难道会产生定性影响？仅仅是带了个nosymm就不一样了，这时候还能够把这种Cs对称性的结果作为下一步分析的结构吗？那么以后如果碰到类似的看起来好像有对称性的结构，应该怎么办呢？比如各种大小的富勒烯这种？
3.对于18碳环，g16中使用小于6-311g*级别的即使构建得不错的def2svp也不行，那么判断能否使用2-zeta基组是不是更多时候还是要更加凭经验判断呢？
4.另外，在当前情况下，虽然都不是D9h对称性，但是wB97XD def2-SVP wB97XD def2-SVP与orca中wB97X-D3结合def2-TZVP的结果都是长短键交替，键角也不是160度，对称性也是Cs，定性一致，所以如果在不提前知道一个分子对称性的情况下，是否有可能出现用大基组的结果再次加剧误导说明小基组的结果的可靠性的可能?

wzkchem5

检查一下是不是默认结构优化收敛限不够严导致的。如果势能面足够平缓，而且结构优化收敛限比较松，那么可能分子还没收敛到足够接近D9h对称性，就满足结构优化收敛要求了。

不想飞的猫头鹰

wzkchem5 发表于 2024-3-18 16:40
检查一下是不是默认结构优化收敛限不够严导致的。如果势能面足够平缓，而且结构优化收敛限比较松，那么可能 ...

嗯嗯，关键词是基于Multiwfn生成的模板修改的，已经按照sob老师的建议加上tightSCF了，但是结果还是如上所述

wzkchem5

不想飞的猫头鹰 发表于 2024-3-18 10:22
嗯嗯，关键词是基于Multiwfn生成的模板修改的，已经按照sob老师的建议加上tightSCF了，但是结果还是如上 ...

我说的是结构优化收敛限，不是SCF收敛限

不想飞的猫头鹰

wzkchem5 发表于 2024-3-18 17:50
我说的是结构优化收敛限，不是SCF收敛限

哦哦，谢啦。我刚才新加了tightopt和verytightopt关键词，在wB97X-D3结合TZVP级别重新计算了一次，关键词分别为! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF tightopt noautostart miniprint nopop和! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF verytightopt noautostart miniprint nopop，分别得到的结构对称性为Cs和C1，仍然没有D9h。为了方便您查看，我把orca的输出文件提供了，可以用ofakeg转换一下。

wzkchem5

不想飞的猫头鹰 发表于 2024-3-18 11:15
哦哦，谢啦。我刚才新加了tightopt和verytightopt关键词，在wB97X-D3结合TZVP级别重新计算了一次，关键词 ...

看起来这个结构离D9h已经非常近了。如果不是因为结构优化收敛限的话，可能是因为格点不够大。用defgrid3甚至refgrid试试？

不想飞的猫头鹰

wzkchem5 发表于 2024-3-18 18:37
看起来这个结构离D9h已经非常近了。如果不是因为结构优化收敛限的话，可能是因为格点不够大。用defgrid3 ...

嗯嗯，昨天看您的回复后加了关键词去试了试，关键词分别为! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF verytightopt noautostart miniprint nopop defgrid3和! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF verytightopt noautostart miniprint nopop refgrid，加了defgrid3后的算完了，加了refgrid的算了很久都没收敛，凌晨四点多停了，刚才又去掉了noautostart尝试续算。加了defgrid3后优化的结构仍然是Cs对称性，输出文件也附上了。

乐平

不想飞的猫头鹰 发表于 2024-3-19 08:29
嗯嗯，昨天看您的回复后加了关键词去试了试，关键词分别为! wB97X-D3 TZVP def2/J RIJCOSX opt freq tigh ...

symm=loose

https://gaussian.com/symmetry/

Tells the program to use looser cutoffs in determining symmetry at the first point. It is designed for use with suboptimal input geometries. Tight says to use the regular criteria at the first point, and it is the default.

不想飞的猫头鹰

乐平 发表于 2024-3-19 09:15
symm=loose

https://gaussian.com/symmetry/

Thanks for your suggestion, but this keyword setting is enabled in Gaussian, which does not exist in ORCA, and symmetry is not enabled in ORCA calculations by default. Using nosymm in Gaussian can disable automatic symmetry. However, just because whether symmetry is enabled or not during the calculation, the resulting structure is qualitatively inconsistent in the judgment of point group symmetry in GaussView, which is my problem.

chrinide

不想飞的猫头鹰 发表于 2024-3-18 18:15
哦哦，谢啦。我刚才新加了tightopt和verytightopt关键词，在wB97X-D3结合TZVP级别重新计算了一次，关键词 ...

这个结构的对称性显然就是D9H，D9H对称的R_cut=0.00239A，这已经是非常小的对称R_cut，如果你非要收敛到D9H对称的R_cut<0.001A，对ORCA这种不支持对称约束的程序而言，直接优化到对称R_cut<0.001A非常难；简单的方法就是用当前收敛的结构生成对称匹配的新坐标，再用ORCA接着优化(可以适当提高相关grid和收敛限)--重复上面的步骤多次(一般2-3次就可以)，应该就可以达到。第二个方法就是用其他支持对称约束程序，比如G16/GAMESS-US/TM用几乎相同的条件优化一个D9H的结构，保存波函数(以及HESS-如果做了频率分析)，转换为ORCA格式再优化/或者直接做解析频率计算，直接就验证Force是否收敛，不出意外基本上一步到位（当然不推荐这么干）

不想飞的猫头鹰

chrinide 发表于 2024-3-19 10:35
这个结构的对称性显然就是D9H，D9H对称的R_cut=0.00239A，这已经是非常小的对称R_cut，如果你非要收敛到D ...

嗯嗯，感谢热心网友。也就是说，在不启用对称性的情况下，想要产生GaussView自动判断成高对称性的结构是非常困难的，这和使用的算法有关吧，那么对于18碳环这种已经证明了有D9h对称性的，还能否用ORCA直接优化+振动分析的无法在GaussView中判断出高对称性的结构进行下一步分析吗？
我按照卢老师博文http://sobereva.com/517《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》操作了一下，将卢老师在http://sobereva.com/515《谈谈18碳环的几何结构和电子结构》提供的在wB97X-D/TZVP级别下的C18.fchk转化成gbw文件后，用同样的结构坐标，关键词! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF miniprint nopop来读取gbw文件，很快就收敛了，但是生成的结构在gview里仍然是判断成了Cs对称性，长短键交替，键角总是差一丁点达到160度，orca的输出文件附带上了，可用ofakeg转换。

wzkchem5

chrinide 发表于 2024-3-19 03:35
这个结构的对称性显然就是D9H，D9H对称的R_cut=0.00239A，这已经是非常小的对称R_cut，如果你非要收敛到D ...

另外ORCA也可以用UseSym关键字限制在某个特定的点群下面优化，而且有%sym模块的symthresh关键字可以调节判断对称性的阈值，把symthresh设松一点有利于程序把这个分子判断成D9h

不想飞的猫头鹰

wzkchem5 发表于 2024-3-19 16:35
另外ORCA也可以用UseSym关键字限制在某个特定的点群下面优化，而且有%sym模块的symthresh关键字可以调节 ...

嗯嗯，感谢您提供的思路，我看了下orca5.0.4的手册，在60/1330页说这个功能不太成熟（THIS IS VERY RUDIMENTARY!），和88/1330页的“ On request, with the UseSym keyword or a %symmetry input block (for which the abbreviation %sym is allowed), the program detects the point group, cleans up the coordinates, orients the molecule and produces symmetry-adapted orbitals in SCF/CASSCF calculations. Note however that the calculation time will not be reduced. While for geometry cleanup the full point group is taken into account (as long as the number of group operations does not exceed 120), only D2h and subgroups are currently supported for electronic-structure calculations. The only correlation module that makes use of this information so far is the MRCI module. Here and in CASSCF calculations, the use of symmetry helps to control the calculation and the interpretation of the results. More symmetry is likely to be implemented in the future, although it is unlikely that the program will ever take advantage of symmetry in a very big way.”说目前不支持太高的对称性。
这种特殊情况也并不好判断，什么时候启用对称性或者用更高质量的格点也有点人为性。如果不按照gview中那么严格的判断对称性的标准的话，就按照卢老师多次调试过的生成的对于orca 5.0以上的关键词，在! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF verytightopt noautostart miniprint nopop关键词下，生成的18个键角为159.86327 159.87491 159.96979 160.07819 160.11089 160.08013 160.08016 160.01877 159.94186 159.86871 159.85146 159.98282 160.01989 160.13435 160.12767 160.07948 160.00316 159.91447 度，而加了defgrid3和verytightopt关键词的18个键角为160.01266 160.01485 160.00878 160.00682 159.99060 159.97376 159.97637 159.98734 160.01363 160.03075 160.02324 160.01787 159.98388 159.97666 159.97801 159.98562 160.00905 160.01013，也只是整体只靠近了均匀的160度一点而已，但是耗时翻了很多倍，之前收敛失败的refgrid重新读取初猜后又收敛失败了。这样的话，之后如果想研究类似的已经证明有对称性的分子，类似18碳环这种，用ORCA进行优化和振动分析是否合适，如果把键长键角列出来后说这已经十分接近D9h了的话是否足够有说服力呢？当前体系在wB97X-D3 TZVP的结果可以用来做分析吗？

wzkchem5

不想飞的猫头鹰 发表于 2024-3-19 10:09
嗯嗯，感谢您提供的思路，我看了下orca5.0.4的手册，在60/1330页说这个功能不太成熟（THIS IS VERY RUDIM ...

我的观点是，可以把分子对称化一下，得到严格D9h对称性的结构，再做一个opt freq，如果没有虚频，那就可以辅助说明真实结构确实是D9h，而不是离D9h特别近的一个非D9h结构。
事实上，如果一个分子的全局最稳定结构几乎是对称的，而且对称结构和全局最稳定结构的能量差远小于kT，那么它和完全对称也没有本质区别了，因为室温下的热平均结构就是对称的。如果对称结构的能量再低一些，可能即使在0K下，也会因为隧穿，导致实验只能看到对称的平均结构。此时去纠结这个分子的平衡结构是完全对称还是几乎对称，可能不是太有意义。

不想飞的猫头鹰

wzkchem5 发表于 2024-3-19 17:15
我的观点是，可以把分子对称化一下，得到严格D9h对称性的结构，再做一个opt freq，如果没有虚频，那就可 ...

嗯嗯，一开始最初使用的结构就是卢老师提供的满足D9h对称性的fchk文件转化出的坐标，然后在orca优化后对称性反而被破坏了。卢老师wb97xd TZVP的fchk文件的能量是-685.222012657429 Hartree，使用这个结构做orca的输入文件，在关键词! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop下，取DFT DISPERSION CORRECTION字段后（整个out文件最后的一个FINAL SINGLE POINT ENERGY字段）得到的能量是-685.256502269139，在关键词! wB97X-D3 def2-QZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop下得到的能量是-685.323267917063，在关键词! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF verytightopt noautostart miniprint nopop defgrid3下得到的能量是-685.256278164891，orca本身计算结果之间的差距不大，跟gaussian的差别应该可以当成D2 D3的校正能与格点设置之类的程序计算方法的差别。
所以说，即使orca的结果没有被gview判断出对称性，也可以放心地提供坐标，以及进一步分析了吧？另外，如果是要在正文中说明结构特征，有必要提一下体系的结构对称性吗？