计算化学公社

标题: 动手尝试18碳环的优化与振动分析出现的问题 [打印本页]

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-18 13:19
标题: 动手尝试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，定性一致，所以如果在不提前知道一个分子对称性的情况下，是否有可能出现用大基组的结果再次加剧误导说明小基组的结果的可靠性的可能?

---

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

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-18 17:22

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

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

---

作者
Author: wzkchem5    时间: 2024-3-18 17:50

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

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

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-18 18:15

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转换一下。

---

作者
Author: wzkchem5    时间: 2024-3-18 18:37

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

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

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-19 08:29

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对称性，输出文件也附上了。

---

作者
Author: 乐平    时间: 2024-3-19 09:15

不想飞的猫头鹰 发表于 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.

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-19 09:53

乐平 发表于 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.

---

作者
Author: chrinide    时间: 2024-3-19 10:35

不想飞的猫头鹰 发表于 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是否收敛，不出意外基本上一步到位（当然不推荐这么干）

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-19 16:17

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转换。

---

作者
Author: wzkchem5    时间: 2024-3-19 16:35

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

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

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-19 17:09

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的结果可以用来做分析吗？

---

作者
Author: wzkchem5    时间: 2024-3-19 17:15

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

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

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-19 18:30

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判断出对称性，也可以放心地提供坐标，以及进一步分析了吧？另外，如果是要在正文中说明结构特征，有必要提一下体系的结构对称性吗？

---

作者
Author: wzkchem5    时间: 2024-3-19 21:57

不想飞的猫头鹰 发表于 2024-3-19 11:30
嗯嗯，一开始最初使用的结构就是卢老师提供的满足D9h对称性的fchk文件转化出的坐标，然后在orca优化后对 ...

可以把orca优化的结构再对称化一下，再用orca优化。
orca和高斯因为格点、泛函/基组细节等不一致，平衡结构会有微小区别。在不用UseSym的情况下，即便初始结构和收敛的结构都严格满足对称性，结构优化过程也可能暂时偏离对称性，接近收敛的时候分子才重新趋于对称，这个是因为ORCA在不打开UseSym时，产生的冗余内坐标不能保证满足对称性。
我觉得如果对称性不是讨论的重点，那么直接提供坐标、分析，也是没问题的。有的情况下，不拿到严格对称的结构，可能讨论的时候措辞会比较别扭（比如不能说所有C-C单键都是1.xxx埃，但如果给个范围，别人又会问你的结果是不是表明C18不是D9h对称性，你又没法确凿地说确实不是），这时可以多花一些精力去验证真实结构是不是对称的，例如按我说的限制D9h对称性优化一个结构，再算频率确认没有虚频，以及看能量是不是比Cs结构的能量低。

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-20 11:33

wzkchem5 发表于 2024-3-19 21:57
可以把orca优化的结构再对称化一下，再用orca优化。
orca和高斯因为格点、泛函/基组细节等不一致，平衡 ...

嗯嗯，这样的话，对于一些有对称性的体系用orca研究还挺麻烦，如果真被问到了可就麻烦了，为什么orca在设计的时候不加入像gaussian那样的支持高对称性的算法呢？orca5.0.4手册中没怎么看懂，说是目前只支持D2h及以下的检测？手册88/1330页有个本应该D3d对称性的例子，加的关键词是%sym SymThresh 1.0e-2 end，这个也没看出来是限制在某个对称性下优化，该如何设置限制在D9h下优化呢？

---

作者
Author: wzkchem5    时间: 2024-3-20 16:58

不想飞的猫头鹰 发表于 2024-3-20 04:33
嗯嗯，这样的话，对于一些有对称性的体系用orca研究还挺麻烦，如果真被问到了可就麻烦了，为什么orca在设 ...

首先，orca并不是不支持高对称性点群，只是目前发行版里支持得不完整。如果分子所属点群不是D2h的子群，那么orca已发行的版本中，暂时不支持在电子结构计算里用这个点群。但是（1）orca仍然支持在结构优化里用这个点群；（2）我们内部的版本已经实现了在电子结构计算里用这个点群，只不过还没有发布而已。所以如果指定在D9h点群下优化C18，那么如果用ORCA 5.0.4，结构优化在D9h群下进行，但SCF和梯度在C2v群下计算；等 ORCA 6.0出来以后，结构优化、SCF和梯度都在D9h群下进行。
其次，至于为什么orca 5.0.4不支持用高阶点群做SCF，是因为高阶点群有高维不可约表示或者复表示。D2h及其子群的所有不可约表示都是一维实表示，Cn、Cnh（n>2）、S2n（n>1）群有一维复表示，其他群都有二维或更高维的不可约表示。处理复表示和高维表示需要额外代码，所以假如一个程序只是近期才开始支持对称性，或者支持对称性不是其开发重点，那么这个程序只支持D2h及其子群而不支持其他群是完全可以理解的。甚至可以说，大多数程序都是先支持D2h及其子群，然后才支持的其他群。
最后，关于如何限定在某个特定群下面优化：6.6.5节（视手册版本而定可能未必是这个章节号）明确说了有PointGroup这个关键词，可以指定点群。此外如果程序在默认的（或你设置的）symthresh下判断出的点群恰好符合你的要求，那么不需要指定PointGroup这个关键词，就可以直接限制在自动判断出来的这个点群下优化

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-20 18:44

wzkchem5 发表于 2024-3-20 16:58
首先，orca并不是不支持高对称性点群，只是目前发行版里支持得不完整。如果分子所属点群不是D2h的子群， ...

嗯嗯，刚才试了一下，关键词如下：
! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop Usesym
%sym
PointGroup "D9h"
end
%maxcore  3000
%pal nprocs  36
end
然而这时候报错[file orca_tools/Tool-Symmetry/pgbase.cpp, line 828]: ERROR (ORCA_SYM): The program cannot continue without a valid atom transfer table.以为是对称性设置太高的问题，但改成D2h也不行。
辛苦您的开发和回答，好期待orca的新版呀，要是能像gaussian那样自动启用对称性加速就好了，不需要的时候像gaussian那样写nosymm。
也就是说，当前如果不能用orca优化出键角严格都是160度的碳环的话，结果的结构还能够用于发表文章不？而如果想要研究高对称性的体系与小分子形成的多聚物间的相互作用，预期此时高对称性分子的对称性会被打破，就像http://sobereva.com/572《探究18碳环独特的分子间相互作用与pi-pi堆积特征》里小分子被吸附进碳环那样，这时候orca的结果就不用担心了吧？

---

作者
Author: wzkchem5    时间: 2024-3-20 19:05

不想飞的猫头鹰 发表于 2024-3-20 11:44
嗯嗯，刚才试了一下，关键词如下：
! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF noautostart min ...

可能因为初始结构在默认symthresh下不符合D9h对称性，可以增大symthresh试试。设置D2h肯定是不行的，因为D2h不是D9h的子群。
如果对称性对你研究的问题来说不重要，那么不严格满足对称性的结果可以用

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-20 19:17

wzkchem5 发表于 2024-3-20 19:05
可能因为初始结构在默认symthresh下不符合D9h对称性，可以增大symthresh试试。设置D2h肯定是不行的，因为 ...

嗯嗯，调回了D9h，刚才设置了SymThresh 10.0e-2关键词，一直增大到10.0e-1都不行。目前对于高对称性的单体只能放弃啦，要是只是研究多体相互作用这样没有什么对称性的应该还是可以的。

---

作者
Author: wzkchem5    时间: 2024-3-20 20:39

不想飞的猫头鹰 发表于 2024-3-20 12:17
嗯嗯，调回了D9h，刚才设置了SymThresh 10.0e-2关键词，一直增大到10.0e-1都不行。目前对于高对称性的单 ...

symthresh不是越大越好，可以看到手册里可能最多也就设到1e-2。首先，如果symthresh设得太粗糙，会导致把点群判断成D18h；其次，ORCA的代码未必是在symthresh大到不合理的量级的时候仍然能正常工作的，所以你去看symthresh设得特别大的计算（比如设为1.0）的输出，会发现程序会报一些warning：

1. Warning (ORCA_SYM): Tolerance set to a very high value (1.0E+00)! Beware of odd results!
   
2. The point group will now be determined using a tolerance of 1.0000e+00.
   
3. Splitting atom subsets according to nuclear charge, mass and basis set.
   
4. Splitting atom subsets according to distance from the molecule's center.
   
5. Identifying relative distance patterns of the atoms.
   
6. Splitting atom subsets according to atoms' relative distance patterns.
   
7. Bring atoms of each subset into input order.
   
8. The molecule is planar.
   
9. The molecule has a center of inversion.
   
10. Analyzing the first atom subset for its symmetry.
    
11. 
    
12. WARNING (ORCA_SYM): Unexpected number of nearest neighbors (17) in function
    
13. TPolytope::CenteredPolygon(TSymmetryFinder &).
    
14. 
    
15. WARNING (ORCA_SYM): None of the atom subsets could be assigned a point group.
    
16. Please take all warnings, especially the ORCA_SYM warnings, very seriously!!!
    
17. 
    
18. WARNING (ORCA_SYM): An inconsistency has occured during symmetry recognition.
    
19. Using C2h as a fallback. If this point group appears to be too low, please try
    
20. increasing the cleaning tolerance using option SymThresh in the %symmetry input
    
21. block. If the problem persists, you may post a message in the ORCA forum.
    
22. 
    
23. Testing point group C2h.
    
24. Success!
    
25. This point group has been found:    C2h
    
26. Largest non-degenerate subgroup:    C2h
    

复制代码

可以看到，首先程序明确说这个symthresh设得太大了，结果可能会比较奇怪，然后报了好几个警告，因为判断点群的时候出现了内部不自洽的问题，所以放弃判断真实点群，但是程序至少知道C2h是分子的子群，所以后续用C2h群计算，虽然未必是用户想要的结果，但是起码能让计算正常进行。
另外我这里测试的结果是，symthresh设成1e-2是可以正确识别点群为D9h的。所以程序应该没问题。做计算一定要看输出文件里的警告，不能只看最终结果

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-20 22:26

wzkchem5 发表于 2024-3-20 20:39
symthresh不是越大越好，可以看到手册里可能最多也就设到1e-2。首先，如果symthresh设得太粗糙，会导致把 ...

嗯嗯，十分感谢！我这里可能还是输入文件有问题，把SymThresh设为1e-2还是出现同样的错误，去看了下输出的out文件发现有警告：
WARNING: Geometry optimizations can change the point-group symmetry!


WARNING: Geometry Optimization
  ===> : Switching off AutoStart
         For restart on a previous wavefunction, please use MOREAD
说是结构优化会改变点群对称性，然后我尝试把它换成单点任务，就终于不报错了，但是奇怪的是，原fchk文件在gview里是D9h对称性，而用orca单点计算后用ofakeg转换出来的fake.out文件在用gview打开的时候一开始会有个警告说任务仍然在跑，可以无视，然后点群对称性判断为D9d，直接把gbw转换出molden再用Multiwfn转成gif，gview判断的点群仍为D9d。
我用来结构优化和振动分析的关键词如下（做单点只是把opt freq去掉了）：
! wB97X-D3 TZVP def2/J RIJCOSX opt freq tightSCF miniprint nopop Usesym
%sym
PointGroup "D9h"
SymThresh 1e-2
end
这个文件里的坐标是直接用卢老师那个fchk的坐标。但是换成普通的orca直接优化振动分析出来的坐标也是这样的错误。输入文件我也附上了，您方便的话可以看下报如上的警告和错误是哪里出问题了吗？

---

作者
Author: wzkchem5    时间: 2024-3-21 01:24

不想飞的猫头鹰 发表于 2024-3-20 15:26
嗯嗯，十分感谢！我这里可能还是输入文件有问题，把SymThresh设为1e-2还是出现同样的错误，去看了下输出 ...

排查了一下，问题出在：如果你要指定点群，那么你的分子必须一开始就处于正确的标准朝向。说明书讲PointGroup关键词的地方写了：

1. If the user specifies a point group using this option, point group
   
2. recognition will be skipped and the user must make sure that
   
3. the molecule is oriented in the coordinate system in agreement
   
4. with the conventions in Section 6.6.3.

复制代码

你的初始结构确实是严格的D9h对称性，但不符合ORCA里面D9h群的标准取向。在5.0.4版里，程序确实有一点做得不尽人意，就是报错退出的时候没有把报错原因说得很明白，我是用内部的orca版本（这里报错信息改进得完善了一些）才发现的问题。
所以由此看来，如果你已知你的分子属于某个点群，而且指定就要用最高可用的点群来做计算，那么最好不在输入文件里指定点群，而是通过适当调整symthresh使得程序能够自动认出来这个点群。要不然的话，程序看不出这个分子哪里像D9h了，你即使硬告诉它这个分子就是D9h，它也不知道这个分子该怎么摆，对称轴、对称面在哪里，而如果程序知道对称轴、对称面在哪里的话，它就看出这个分子是D9h了，不需要你告诉它。输入文件里指定点群的做法主要是为了让你可以指定用一个子群做计算，例如真实对称性D9h，但你只需要保持C2v对称性，那就可以指定C2v。

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-21 17:15

wzkchem5 发表于 2024-3-21 01:24
排查了一下，问题出在：如果你要指定点群，那么你的分子必须一开始就处于正确的标准朝向。说明书讲PointG ...

嗯嗯，按照您说的，用Multiwfn的几何操作功能把碳环放到了垂直于或者平行于XY平面的方向，并且去掉了point group那行关键词，单点任务出来的结构终于保持了原有的D9h对称性。不过opt freq任务仍然会生成Cs对称性，看来目前如果怀疑一个体系有高对称性，就只能用高斯优化后，orca做单点计算了，感觉如果分子本身结构复杂不好考虑应该怎么摆放，单点计算那一点误差应该无伤大雅。
这个过程太曲折了，哈哈，orca能够把对称性的功能完善一下就太好了。也感谢您的不懈耐心回答和对orca的坚持完善！

---

作者
Author: wzkchem5    时间: 2024-3-21 18:12

不想飞的猫头鹰 发表于 2024-3-21 10:15
嗯嗯，按照您说的，用Multiwfn的几何操作功能把碳环放到了垂直于或者平行于XY平面的方向，并且去掉了poin ...

刚发现问题没有我想象中的那么简单，因此编辑了一下帖子。
我发现C18在wB97X-D3/TZVP级别下的波函数其实是不满足D9h对称性的。关键词：
! wB97X-D3 TZVP def2/J RIJCOSX extremeSCF Usesym defgrid3 autostart
其中坐标是由ORCA先做一个不指定pointgroup的计算，得到的旋转到标准取向的D9h结构。格点加到defgrid3，以排除格点的影响；收敛限加到extremeSCF，以排除没有完全收敛的影响。
可以看到收敛的波函数破坏了D9h对称性，碳的Mulliken电荷略微偏离0：

1.    0 C :   -0.009357
   
2.    1 C :   -0.019616
   
3.    2 C :    0.038474
   
4.    3 C :    0.048717
   
5.    4 C :   -0.036005
   
6.    5 C :   -0.022217
   
7.    6 C :   -0.022217
   
8.    7 C :   -0.036005
   
9.    8 C :    0.048718
   
10.    9 C :    0.038474
    
11.   10 C :   -0.019616
    
12.   11 C :   -0.009357
    
13.   12 C :    0.002223
    
14.   13 C :   -0.028843
    
15.   14 C :    0.026623
    
16.   15 C :    0.026623
    
17.   16 C :   -0.028844
    
18.   17 C :    0.002224

复制代码

因为TRAH内置了通过随机扰动打破对称性的功能，所以可以找到能量更低的、破坏了点群对称性的波函数。那么既然波函数的对称性低于D9h，梯度的对称性必然低于D9h，所以结构优化必然会破坏D9h对称性，不破坏反而是错的。
可以检查一下你的高斯计算的波函数是否严格满足D9h对称性，如果严格满足，检查波函数是否有RKS instability。搞不好高斯收敛到的波函数不是稳定波函数。
总之结论就是，结构优化即使指定了点群，也未必不应该破坏对称性，因为波函数的对称性可能低于分子结构的对称性。

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-21 22:06

wzkchem5 发表于 2024-3-21 18:12
刚发现问题没有我想象中的那么简单，因此编辑了一下帖子。
我发现C18在wB97X-D3/TZVP级别下的波函数其实 ...

嗯嗯，fchk文件就是卢老师在http://sobereva.com/515《谈谈18碳环的几何结构和电子结构》里提供的wB97XD.fchk文件，应该是稳定的，用记事本打开也没有instability字样。
至于即使只用orca做单点的情况，如您所述，其实orca的单点计算方法由于会故意打破波函数对称性，所以这时候的坐标也会轻微变动，导致单点出来的结构在gview里仍然只有Cs对称性了？也就是说，orca目前的算法就不太适用于优化明确具有高对称性的体系？

---

作者
Author: wzkchem5    时间: 2024-3-22 01:03

不想飞的猫头鹰 发表于 2024-3-21 15:06
嗯嗯，fchk文件就是卢老师在http://sobereva.com/515《谈谈18碳环的几何结构和电子结构》里提供的w ...

哦，sob老师用的是wB97XD，咱们的orca计算用的是wB97X-D3，这两者是不同的泛函，即便去掉色散校正部分以后也不完全相同。虽然orca不支持wB97XD的色散校正部分，但是支持wB97XD除去色散校正以外的部分（需要用输入libxc泛函的语法来输入）。我试了一下wB97XD（不带色散校正部分），用ORCA 5.0.4跑的：

! TZVP def2/J RIJCOSX tightSCF defgrid3 opt bohrs
%method
method dft
functional hyb_gga_xc_wb97x_d
end
%sym
SymThresh 1e-2
end
%pal nprocs  8 end
* xyz 0 1
C     2.331446759384    6.577411361089    0.000000000000
C     0.058763860883    6.978146674690    0.000000000000
C    -2.441886692266    6.537214512846    0.000000000000
C    -4.440450491980    5.383343165088    0.000000000000
C    -6.072634222057    3.438182340995    0.000000000000
C    -6.861928709536    1.269613559347    0.000000000000
C    -6.861928709536   -1.269613559347    0.000000000000
C    -6.072634222057   -3.438182340995    0.000000000000
C    -4.440450491980   -5.383343165088    0.000000000000
C    -2.441886692266   -6.537214512846    0.000000000000
C     0.058763860883   -6.978146674690    0.000000000000
C     2.331446759384   -6.577411361089    0.000000000000
C     4.530481950152   -5.307797801742    0.000000000000
C     6.013870361173   -3.539964333694    0.000000000000
C     6.882337184246   -1.153871347758    0.000000000000
C     6.882337184246    1.153871347758    0.000000000000
C     6.013870361173    3.539964333694    0.000000000000
C     4.530481950152    5.307797801742    0.000000000000
*

结果发现是可以得到满足D9h对称性的波函数的，而且一步就收敛了。收敛的能量-685.215436111572比fchk文件里略高（-6.852220126574286E+02），和我这个计算没加色散校正、而sob老师加了也是吻合的。充分说明orca不仅在波函数该对称的时候能保持对称性，而且计算结果和高斯足够接近（因为这个结构是sob老师用高斯优化过的，虽然受orca程序限制没办法加wB97XD的D2色散校正，但这个体系不太大且不太柔性，色散校正不重要）。
另外，如果已知一个体系肯定有高对称性，但是orca可以收敛到低对称性的波函数，而且能量比高对称性波函数更低，那么说明用的泛函不适合（该泛函收敛得到的最稳定的波函数的对称性与实际不符），而不是orca有问题。因为如果当前泛函的最稳定波函数是具有高对称性的，那么orca即使暂时破坏了对称性（不管是有意还是无意），最后也能收敛回来。这个和计算开壳层单重态是一样的道理，如果体系唯一的稳定波函数是闭壳层单重态，你做一个broken symmetry计算，虽然一开始会破坏自旋对称性，但是最后这个开壳层单重态会收敛回闭壳层单重态。

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-22 10:46

wzkchem5 发表于 2024-3-22 01:03
哦，sob老师用的是wB97XD，咱们的orca计算用的是wB97X-D3，这两者是不同的泛函，即便去掉色散校正部分以 ...

嗯嗯，我在帖子一开始把关键词和计算级别写的比较清楚，最初就是想试一下本以为更好的wb97x-d3，但是目前看来，wb97x-d3并没有产生正确的结构。
在卢老师博文http://sobereva.com/272《简谈量子化学计算中DFT泛函的选择》说“wB97XD提出后，其中作者之一蔡政达又提出了wB97X-D3，把色散校正从DFT-D2改成了DFT-D3，也改了些其它参数。wB97X-D3的精度比wB97XD整体有所改善，包括热化学、弱相互作用等方面，只要程序能支持wB97X-D3（ORCA支持，但Gaussian不支持），那么前面说的wB97XD都建议改为wB97X-D3。”
看来，应该是“也改了些其他参数”这部分对计算结果产生了定性影响。也就是说，并不是任何情况下，都应该用看似更好的wB97X-D3？
我将您的关键词和坐标复制，尝试了下优化任务，最后的结构取出来用ofakeg处理，在gview里发现current point group仍然是Cs，但是aproximaet higher-order point group这里判断出了D9h。如果把opt关键词去掉，做单点计算，这次终于在current point group里也是D9h了，而在aproximaet higher-order point group没有字样，把右边可选的tolerance从默认的0.01改为0.1，会被判断为D9d。
同样的初始坐标下，换了个目前表现最好的非双杂化泛函wb97m-v, opt的结构current point group里是C1，aproximaet higher-order point group没有字样，tolerance改为0.1后是D9d，而去掉opt改成单点计算，结果是current point group里是D9h，tolerance改为0.1时仍然不显示，0.3后才是D9d，感觉这里的gview判断标准有点迷，approximate higer-order point groups判断的不一定准，current point group的也可能有问题。
还有如前文所述，尝试用wb97x-d4（orca手册上说这个是wb97m-v修改而来）优化的结果也是没有对称性，把已经修改好坐标使碳环处于XY平面的上的坐标作为初始坐标以满足对称性检测的要求，直接就SCF不收敛了，opt或者单点都是，即使加上defgrid3也不行。又试了试M062X，也是不收敛。
也就是说，这个问题还涉及到泛函本身的问题。不知之后的orca未来可否加上与gaussian相同定义的D2版本的wb97xd，以及对于不依赖于初始摆位的对称性的启用，来和gaussian对比呢？

---

作者
Author: wzkchem5    时间: 2024-3-22 17:23

不想飞的猫头鹰 发表于 2024-3-22 03:46
嗯嗯，我在帖子一开始把关键词和计算级别写的比较清楚，最初就是想试一下本以为更好的wb97x-d3，但是目前 ...

wB97X-D3的色散校正部分，确实本来就该比wB97XD好。但是其他部分未必，尤其是wB97X-D3的HF exchange比例和wB97XD不同。众所周知，有的体系适合HF exchange比例小，有的体系适合HF exchange比例大，所以如果改了HF exchange比例，必然会导致对有的体系精度变差，最多只能希望变差的情况比变好的情况少。
ORCA的色散校正只实现了D3和D4，没有计划实现较老的、精度较差的D2。固然这会使得难以重复别人做的wB97XD计算，但是正如我所说，wB97XD即使比wB97X-D3好，也一般不是因为它用了D2，而是它的HF exchange比例碰巧适合这个体系。那完全可以找另外一个HF exchange差不多的泛函去做计算。如果只是为了和高斯对比，应该也可以在高斯里把wB97XD的色散校正关掉。
SCF不收敛可能是因为波函数恰好在破缺对称性和不破缺对称性的边缘，此时波函数对于外界扰动极其敏感。可以想象一下你优化f(x) = x^4 + kx^2这个函数，当k接近0的时候，极小值点处的二阶导也接近0，函数曲线极其平缓，且非谐性严重，导致一般的优化算法收敛变慢。而k=0恰好是这个函数极小值点偏离对称性的临界点

---

作者
Author: 不想飞的猫头鹰    时间: 2024-3-25 20:09

wzkchem5 发表于 2024-3-22 17:23
wB97X-D3的色散校正部分，确实本来就该比wB97XD好。但是其他部分未必，尤其是wB97X-D3的HF exchange比例 ...

不好意思啦，前两天出去了。非常感谢您提供的视角，泛函的选取和特点果然得多测试和理解泛函提出的形式与目的，看来对于当前的体系D2形式的wb97xd对于碳环的描述正好比较沾光。再次感谢您的耐心回答，这对未来类似体系泛函的考虑与选取十分有帮助。

---

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

Powered by Discuz! X3.3