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

wzkchem5

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

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

不想飞的猫头鹰

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下优化呢？

wzkchem5

不想飞的猫头鹰 发表于 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这个关键词，就可以直接限制在自动判断出来的这个点群下优化

不想飞的猫头鹰

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的结果就不用担心了吧？

wzkchem5

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

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

不想飞的猫头鹰

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

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

wzkchem5

不想飞的猫头鹰 发表于 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的。所以程序应该没问题。做计算一定要看输出文件里的警告，不能只看最终结果

不想飞的猫头鹰

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直接优化振动分析出来的坐标也是这样的错误。输入文件我也附上了，您方便的话可以看下报如上的警告和错误是哪里出问题了吗？

wzkchem5

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

不想飞的猫头鹰

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

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

wzkchem5

不想飞的猫头鹰 发表于 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。搞不好高斯收敛到的波函数不是稳定波函数。
总之结论就是，结构优化即使指定了点群，也未必不应该破坏对称性，因为波函数的对称性可能低于分子结构的对称性。

不想飞的猫头鹰

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

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

wzkchem5

不想飞的猫头鹰 发表于 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计算，虽然一开始会破坏自旋对称性，但是最后这个开壳层单重态会收敛回闭壳层单重态。

不想飞的猫头鹰

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对比呢？

wzkchem5

不想飞的猫头鹰 发表于 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恰好是这个函数极小值点偏离对称性的临界点