计算化学公社

标题: 求助：自由能计算的时候，怎样禁止ORCA考虑对称性 [打印本页]

作者
Author: dnlx 时间: 2022-4-25 11:08
标题: 求助：自由能计算的时候，怎样禁止ORCA考虑对称性
本来以为ORCA的计算从来不考虑对称性，但今天发现不是这样。我在计算水溶液中乙烷的自由能G-E(el)，三种碳同位素取代情况：没有13C，一个13C，和两个13C。此前用B3LYP-D4/dev2-SVP优化的时候，模拟的结果没有包括分子对称性的因素，G-E(el)的次序是：

没有取代的乙烷 (25.02 kcal/mol) > 一个13C取代的乙烷 (24.85 kcal/mol) > 两个个13C取代的乙烷 (24.69 kcal/mol)

现在用B3LYP-D4/dev2-TZVP优化（关键词 !B3LYP D4 def2-TZVP def2/J RIJCOSX Opt VeryTightOPT strongSCF NumFreq nopop，因为在水溶液中所以用NumFreq），发现模拟结果包含了对称性（从取代位置以及反应速度上推断），三者的G-E(el)是

24.45 kcal/mol < 24.86 kcal/mol > 24.12 kcal/mol

现在的问题是：
一、为什么ORCA （5.0.1）的计算结果，有时候包括了对称性，有时候不包括？
二、怎样得到不包含对称性的计算结果？

请高手指点！

作者
Author: fatpig 时间: 2022-4-25 12:04
看看你的output file里有没有“rotational symmetry number”(Sn)，当Sn=1时就是C1也就是无对称性的结果。
如果没有的话，应该可以在一开始就指定用无对称性进行计算。
如果你得到了有对称性的结果，也可以找出这个对称性的symmetry number(Sn)，然后用公式换算成Sn=1的结果。

For a nonlinear molecule the correct rotational entropy is:
S(rot) = R*(ln(qrot/sn)+1.5)
R = 8.31441 J/mol/K = 1.987191683e-3 kcal/mol/K
qrot = 4922043.1563766
sn is the rotational symmetry number. We have assumed 3 here
   if it is different for your molecule then you should correct
   the printed rotational entropy by manually evaluating the equation
   as given above

For convenience we print out the resulting values for sn=1 - 12:
sn= 1  qrot/sn= 4922043.1564 T*S(rot)= 10.02 kcal/mol T*S(tot)= 34.72 kcal/mol
sn= 2  qrot/sn= 2461021.5782 T*S(rot)=    9.61 kcal/mol T*S(tot)= 34.31 kcal/mol
sn= 3  qrot/sn= 1640681.0521 T*S(rot)=    9.37 kcal/mol T*S(tot)= 34.07 kcal/mol
sn= 4  qrot/sn= 1230510.7891 T*S(rot)=    9.20 kcal/mol T*S(tot)= 33.90 kcal/mol
sn= 5  qrot/sn=  984408.6313 T*S(rot)=    9.06 kcal/mol T*S(tot)= 33.77 kcal/mol
sn= 6  qrot/sn=  820340.5261 T*S(rot)=    8.96 kcal/mol T*S(tot)= 33.66 kcal/mol
sn= 7  qrot/sn=  703149.0223 T*S(rot)=    8.87 kcal/mol T*S(tot)= 33.57 kcal/mol
sn= 8  qrot/sn=  615255.3945 T*S(rot)=    8.79 kcal/mol T*S(tot)= 33.49 kcal/mol
sn= 9  qrot/sn=  546893.6840 T*S(rot)=    8.72 kcal/mol T*S(tot)= 33.42 kcal/mol

作者
Author: wzkchem5 时间: 2022-4-25 15:30
可能是你一开始画的分子不是严格对称的，第一个计算的结构优化收敛限比较粗糙，没有把不对称的结构优化成对称的结构，但是第二个计算的收敛限比较严，把结构优化成完全对称的了

作者
Author: dnlx 时间: 2022-4-25 20:48
本帖最后由 dnlx 于 2022-4-25 20:54 编辑

wzkchem5 发表于 2022-4-25 15:30
可能是你一开始画的分子不是严格对称的，第一个计算的结构优化收敛限比较粗糙，没有把不对称的结构优化成对 ...

第一个确实只用了默认的收敛限。但是麻烦在于第二个，其输出文件给的是“Point Group: Ci, Symmetry Number: 1 ”，然而替代了一个同位素之后，识别出的对称性反而增加了“Point Group: C2, Symmetry Number: 2 ”

作者
Author: wzkchem5 时间: 2022-4-25 20:57

dnlx 发表于 2022-4-25 13:48
第一个确实只用了默认的收敛限。但是麻烦在于第二个，其输出文件给的是“Point Group: Ci, Symmetry Num ...

可能是替代同位素的那个计算碰巧收敛得更严了一些？

作者
Author: dnlx 时间: 2022-4-25 21:06

wzkchem5 发表于 2022-4-25 20:57
可能是替代同位素的那个计算碰巧收敛得更严了一些？

同位素替换之后的计算只是热力学计算，不明白为什么会影响对称性的识别。下面是相关的输出文件。

(, 下载次数 Times of downloads: 10) (, 下载次数 Times of downloads: 5) (, 下载次数 Times of downloads: 4) (, 下载次数 Times of downloads: 1)

作者
Author: wzkchem5 时间: 2022-4-25 22:05

dnlx 发表于 2022-4-25 14:06
同位素替换之后的计算只是热力学计算，不明白为什么会影响对称性的识别。下面是相关的输出文件。

看起来确实是程序的问题，我跟负责这块程序的人反映一下

作者
Author: wzkchem5 时间: 2022-4-26 21:25

dnlx 发表于 2022-4-25 14:06
同位素替换之后的计算只是热力学计算，不明白为什么会影响对称性的识别。下面是相关的输出文件。

问到了，确实是程序的bug，对于某些同位素取代的分子，点群会识别错。下一版orca会修复。
在此之前如果做同位素取代的、高对称性分子的频率计算，应当检查rotational symmetry number是否正确，如果不正确的话按2楼的方法手动计算正确的转动熵。
此外针对这个特定的分子而言，默认的判断对称性的阈值还是太严了，没有办法判断出分子的三重轴。所以还要检查点群有没有识别对，如果识别出来的点群比实际点群的对称性低，需要设一个比较松的阈值，比如%sym symthresh 1e-3 end，然后重新跑（这一点和bug无关，是本来用户就应该检查的）。

作者
Author: Arana 时间: 2023-2-5 12:22
请问下ORCA同位素设置具体是怎么做的呢？查阅手册是这么写的isotopes or nuclear charges are specified with the statements “M = . . . ” and “Z =
. . . ”, respectively, after the atomic coordinate definition.因为所研究的是氘代物质，在原子坐标后添加M=2作结构优化，但是输出文件仍显示原子为H，之前sobl老师讲同位素不影响结构，但输出文件为什么原子还是显示H呢，是哪里有问题嘛

作者
Author: wzkchem5 时间: 2023-2-5 16:25

Arana 发表于 2023-2-5 05:22
请问下ORCA同位素设置具体是怎么做的呢？查阅手册是这么写的isotopes or nuclear charges are specified wi ...

因为D也叫2H（2上标），所以显示成H没有问题。
其他同位素都不会像氘、氚这样有单独的元素符号，所以程序肯定要默认在同位素取代以后元素符号不变。那么问题就是要不要考虑氘、氚的特殊情况，如果显示成D、T，好处是用户看起来没有歧义；如果显示成H，好处是其他的可视化程序、分析计算结果的脚本等可以确保正常读取这个文件（如果显示成D，搞不好有的考虑不周的可视化软件就打不开这个输出文件了，告诉你D这个元素不存在，输出文件有误）。ORCA团队经过权衡以后，认为后者更重要，所以D、T一律显示为H。

作者
Author: Arana 时间: 2023-2-6 08:33

wzkchem5 发表于 2023-2-5 16:25
因为D也叫2H（2上标），所以显示成H没有问题。
其他同位素都不会像氘、氚这样有单独的元素符号，所以程 ...

好的，感谢老师答复。但是在我使用multiwfn帮助生成ORCA输入文件时，读氘代化合物结构优化完的xyz文件，分子式会显示H，而且分子量也是氢的分子量，这个是不是有问题。不好意思，占了11#

作者
Author: sobereva 时间: 2023-2-6 08:42

Arana 发表于 2023-2-6 08:33
好的，感谢老师答复。但是在我使用multiwfn帮助生成ORCA输入文件时，读氘代化合物结构优化完的xyz文件， ...

Multiwfn的这个功能不考虑同位素
自己修改生成的ORCA输入文件，添加非标准同位素设置就完了

作者
Author: Arana 时间: 2023-2-6 08:46

sobereva 发表于 2023-2-6 08:42
Multiwfn的这个功能不考虑同位素
自己修改生成的ORCA输入文件，添加非标准同位素设置就完了

好的，感谢sob老师

作者
Author: wzkchem5 时间: 2023-2-6 16:11

Arana 发表于 2023-2-6 01:33
好的，感谢老师答复。但是在我使用multiwfn帮助生成ORCA输入文件时，读氘代化合物结构优化完的xyz文件， ...

一样的原因，因为显示结构不需要知道用的是什么同位素，而且很多可视化程序也不支持显示同位素，所以xyz文件里根本就没有同位素信息

作者
Author: Arana 时间: 2023-2-8 20:58

wzkchem5 发表于 2023-2-6 16:11
一样的原因，因为显示结构不需要知道用的是什么同位素，而且很多可视化程序也不支持显示同位素，所以xyz ...

好的，谢谢老师

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