计算化学公社
标题: 谈谈Gaussian中的对称性与nosymm关键词的使用 [打印本页]
作者Author: sobereva 时间: 2015-6-27 19:17
标题: 谈谈Gaussian中的对称性与nosymm关键词的使用
谈谈Gaussian中的对称性与nosymm关键词的使用
On the symmetry in Gaussian program and the use of nosymm keyword
文/Sobereva @北京科音 First release: 2015-Jun-27 Last update: 2021-Nov-21
经常有人问Gaussian中的nosymm关键词怎么用、什么时候用，网上也有很多帖子、资料都涉及nosymm，但是很多都明显错误，误人子弟。还有很多人乱用nosymm，不仅没什么意义还白浪费了大量时间。故写此文澄清一下nosymm的用法，以正视听。
1 关于利用对称性
为了更好地理解nosymm，首先说说对称性问题。
没对称性的分子所属点群就是C1，有对称性的分子则属于比C1更高阶的点群。量化计算中电子积分的计算是很耗时的部分，对于有对称性的分子，要算的电子积分中会有很多是等价的，利用对称性的话就可以避免计算大量重复的积分，从而大大节约时间。一些计算在其它部分也能利用对称性显著节约时间，如积分变换、导数计算等，而且不光能节约时间，很多时候还能显著降低对内存/硬盘的占用。对称性越高的体系，利用对称性可以带来越多的加速，例如用XEON E5-2650v2 16核+G09 D.01在CCSD(T)/def2-TZVP级别下计算面对面构型的苯二聚体单点能，在实际的D6h点群下（这属于很高阶的点群）计算只用一个半小时，而不启用对称性的话则需要多达两天半才能算完。
启动对称性还有个好处是Gaussian可以显示出轨道、电子态、振动模式的不可约表示。如果不启用对称性，即总是当成C1对称性来算，即便体系本身有对称性，不可约表示也会都显示为A，此时轨道的不可约表示就只能基于等值面图判断、电子态的不可约表示需要自己再对占据轨道的不可约表示做直积得到、振动模式的不可约表示也只能基于正则振动坐标的箭头图来判断了，这就麻烦多了。另外，允许程序利用对称性，且正确判断出了实际点群的时候，还可以确保不同不可约表示表示的分子轨道在SCF过程中不发生混合，这对于考察一些问题以及SCF收敛有时会有帮助。
由于利用对称性时有上述优点，所以Gaussian默认是启动对称性的。计算一开始时会自动根据一定算法和阈值判断输入文件里分子的点群。比如苯分子，如果输入文件里的结构完全或几乎满足D6h对称性，那么Gaussian就会把它判断为D6h对称性，并按照D6h的情况来做后续计算。如果输入文件里的结构偏离D6h略明显，但接近其它更低阶的点群，就会被判断为更低阶的点群，比如C2v，也可能完全判断不出对称性，即判断为C1。
输入文件里的坐标称为输入朝向（Input orientation）。Gaussian为了判断和利用对称性，会先把输入文件里的坐标调整到标准朝向（Standard orientation）下，然后再进行点群的判断及后续的各种计算。这个调整过程会把体系进行平移，使原子核电荷中心处于笛卡尔坐标原点，并且根据一定规则对体系进行旋转，具体规则这里不细谈，感兴趣的读者可以看手册里Standard Orientation Conventions这一节的介绍。例如，在计算苯分子的时候，如果你在gview里把分子摆成歪斜着并且中心偏离原点，若在计算过后用gview打开输出文件，就会看到分子已经处在了XY平面，并且苯分子中心被挪到了笛卡尔坐标原点，这就是苯分子在标准朝向下的坐标。值得一提的是，哪怕体系没有对称性，即C1点群，Gaussian也会照样把体系弄到标准朝向下导致体系被平移和旋转。下面是一个示意图
(, 下载次数 Times of downloads: 30)
将体系弄到标准朝向会带来一些麻烦的问题。有些时候，我们必须要求体系的位置和朝向是和输入文件里一致的。比如，我们要沿着某个键的键轴加用field关键词外电场，我们设好了外场矢量，但是如果体系朝向被Gaussian旋转到标准朝向了，那么外电场的方向就不合我们的意愿了（除非输入文件里的坐标已经是标准朝向的坐标），显然结果无意义。再比如，我们按照《使用Multiwfn作电子密度差图》（http://sobereva.com/113）用Multiwfn作二聚体和两个单体之间的密度差图来考察单体形成复合物后电子是怎么转移的，这就要求复合物里面的坐标和单体中的坐标完全一致。如果在计算中允许Gaussian自动把体系平移旋转，那么最后得到的单体及复合物的.wfn/wfx或.fch文件中，单体的原子坐标和它在复合物里的坐标就不对应了，显然密度差图也就毫无意义了。Multiwfn里作CDA分析（《使用Multiwfn做电荷分解分析(CDA)、绘制轨道相互作用图》http://sobereva.com/166）和ETS-NOCV分析（《使用Multiwfn通过ETS-NOCV方法深入分析片段间的轨道相互作用》http://sobereva.com/609）也同样要求整体和片段中的原子坐标必须对应。
2 nosymm的用处
理解了上面关于对称性的讨论，就很容易理解nosymm的带来的好处和坏处了。nosymm的含义是：让Gaussian在计算中完全不利用对称性。
nosymm的用处1：用了nosymm后Gaussian就不会再试图去判断和利用对称性，因此也就不会把体系搞到标准朝向下，这就确保了实际计算过程中的坐标和输入文件里的坐标一致，从而可以避免前面说的诸如外加电场、作密度差图、CDA、ETS-NOCV分析时遇到的问题，所以做这些任务的时候总要加上nosymm。
nosymm的用处2：做几何优化、柔性扫描的时候，Gaussian会把每一步的坐标都调整到标准朝向，这往往导致笛卡尔坐标变化不连续（但内坐标变化是连续的），导致在GaussView里观看优化、扫描轨迹过程中体系朝向会发生剧烈跳变，从而难以考察结构变化过程。解决方法有二，其一就是用nosymm，这样优化的每一步的结构都不会被弄到标准朝向下，而只会输出Input orientation坐标，gview也只会读这些坐标，优化/扫描过程的轨迹看起来就连续了。其二就是不让gview去读取每一步的standard orientation坐标（默认情况）而改为读取input orientation坐标，方法见此文《观看Gaussian优化轨迹时避免结构跳变的方法》（http://sobereva.com/289）。
nosymm的用处3：nosymm还有个情况也应该加，也就是做对称破缺计算的时候（详见《谈谈片段组合波函数与自旋极化单重态》http://sobereva.com/82），如果体系本身是有对称性的，那么加上nosymm往往才能得到对称破缺态。比如将乙烯扭转到90度让pi键被破坏，此时是双自由基状态，如果直接用UB3LYP/6-31G* guess=mix可能因为波函数对称性的约束还是得到闭壳层波函数，而加了nosymm就可以得到双自由基状态了。
nosymm只有以上三个用处，如果你的情况不属于这三个，就别用nosymm！！！对于有对称性的体系，用nosymm带来的坏处是很明显的，首先是不能利用对称性来加速计算，其次是没法显示不可约表示。有些人不懂nosymm有什么用就总是盲目地写上nosymm，往往没带来任何益处，当体系有对称性时还因此白浪费了大量计算时间。我总是向初学者强调：如果不清楚某些关键词是干什么的，就一律不要写！！！没看过此文的话建议看看：《常见的多余的和被滥用的Gaussian关键词》（http://sobereva.com/331）。
3 关于几何优化中维持/破坏对称性
经常有人问在Gaussian的几何优化中怎么维持对称性、怎么破坏对称性，这里明确说一下。首先弄清楚两个事实：
（1）在默认的允许利用对称性时，即便优化一开始把体系判断成了高对称性，优化过程中也可能对称性发生降低甚至变为完全没有对称性。原因要么是因为当前计算级别下低对称性的能量更低，有可能由于一些数值层面的问题破坏了对称性（用弥散函数时容易出现此问题，DFT格点积分精度不够高时也可能会出现此问题，尤其是对M06-2X这种对积分格点精度要求高的泛函来说。参考《密度泛函计算中的格点积分方法》http://sobereva.com/69）。
（2）在使用nosymm禁止使用对称性时，一开始的高对称性结构也依然可能继续保持下去。例如，计算D3h对称性的NH3（原子都处在同一平面，呈正三角形），即便写了nosymm，最后优化出来还是D3h结构，这是因为体系中原子都处在一个平面上，优化过程中受到向平面外的力的分量为0，因此会平面结构会一直保持。
所以，对称性是否会在优化中改变，或者说对称性是否会被限制在初始对称性下，和是否使用nosymm没必然的关系。
如果你想在优化过程中100%保证对称性肯定不发生改变，唯一做法是用内坐标描述结构，而且根据对称性来写内坐标，使得对称等价的坐标共享相同的变量（比如苯，让六个C-C键都共享一个键长变量），且优化的时候用opt=z-matrix让Gaussian在内坐标下优化。当然，这么按照对称性写内坐标对于原子数较多的体系是很麻烦的事情，往往还要借助虚原子。PS：虽然Gaussian里有IOp(2/16=2)可以让优化过程中对称性发生改变时沿用旧对称性，但实际发现不怎么奏效。
对于有对称性的体系的计算，一般的考虑方式是：先让初始结构在gview里做对称化成为你期望的点群，然后用Gaussian优化，查看输出文件确保任务最开始时判断对了实际点群（如果没判断对，用symm=loose）。如果优化过程中点群发生了不希望的降低，且你用了弥散函数，把弥散函数去掉重新优化看看；对于DFT，可以进一步提升积分格点精度再试（如G09的用户可以加int=ultrafine试试，G16的用户可以加int=superfine试试）。如果重新优化后对称性还是出现了下降，换成其它你觉得对此体系也靠谱的泛函，也可以尝试不依赖于积分格点的MP2（前提是适合当前体系）或CCSD（量力而行）。如果对称性还是下降了，那么有极大的可能性高对称性结构确实是不稳定的，最后收敛到的低对称性的结构才是真正的极小点，因此此时若非要保持高对称性则结果没有物理意义。如果最后虽然判断出的对称性下降了，但肉眼观看结构时发现和高对称性的结构几乎看不出什么区别，那么你就当这个结构是高对称性的结构即可，写文章时也可以说它是有高对称性的，对称性的下降可能只是一些很trivial的数值因素引起结构的十分细微的变化，无需在意（之后你甚至可以再做个对称化变成高对称性后直接算后续的单点之类任务）。
如果你刻意允许优化过程中对称性发生变化从而避免被限制在初始的对称性（因为初始的对称性可能不一定是实际极小点结构所具有的），一般做法是人为破坏初始结构中的对称性。比如一个平面体系，你怀疑实际结构不是平面的，那你就随便稍微扭曲其中的一个二面角，或者利用《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（http://sobereva.com/610）里介绍的功能产生对所有原子随机位移后的结构，这样初始结构就没有对称性了，优化时也不可能被强行维持在有对称性的结构上。如果优化后肉眼上看结构又自发变得有对称性了，那就是说明此体系在当前级别下极小点结构就是有对称性（如果你是讲究的人，建议对称化后再重新做优化，结构精度会稍微更高，还能同时给出不可约表示）。另外，你也可以先对有对称性的结构直接做几何优化+振动分析，如果发现优化后的结构仍有对称性但是有虚频，且虚频的振动模式是会导致对称性方式破坏的，就按照虚频模式稍微调节一下结构然后继续优化，通常得到的就是无虚频且对称性比之前更低的实际极小点结构了。
4 其它问题
有人问symm=loose、veryloose是干嘛的。前面提到了，Gaussian会通过一定算法和阈值来根据输入文件里的坐标判断体系的点群。默认的判断标准比较严，即坐标稍微有些偏离点群对称性，就认不出来相应的点群。这有时候会带来些麻烦。比如说，在gview里，选择自带的库里的C60，将它对称化后会看到被gview识别成了Ih点群，但是保存成Gaussian输入文件，然后在Gaussian计算时就会发现被Gaussian判断成了C1点群（也可能判断成别的，和版本有关）。这是因为输入文件里的笛卡尔坐标数值精度有限，我们需要用稍微宽松一点的判断标准，即symm=loose关键词，才能让Gaussian判断出它的实际点群Ih，此时计算速度和C1时完全不在一个数量级。symm=veryloose代表使用更宽松的判断标准。
另一个流行的量子化学程序ORCA不会利用对称性加速计算，因此也不会像Gaussian那样把体系弄到标准朝向下，因此没有类似nosymm这种用途的关键词。也有其它一些程序有类似于Gaussian的做法，比如PSI4程序也会自动把体系搞到标准朝向下，如果想实现nosymm的效果，需要写nocom要求不允许平移，以及noreorient不允许旋转朝向，例如：
molecule MOL {
nocom
noreorient
C 0.00000000 -0.74034185 0.01363083
H 0.87837906 -1.12996571 -0.50403292
...
}
作者Author: jiangning198511 时间: 2015-6-27 20:50
好贴，但是有时优化完后，体系的轨道会出现一些问号，导致最终的态的对称性无法判断。
作者Author: brothers 时间: 2015-6-27 23:59
好帖~好在没用错
作者Author: sobereva 时间: 2015-6-28 05:41
即便判断对了点群，也可能有问号出现，这是高斯自己的问题（不同版本情况可能不一样，比如g03没判断出来时可能g09就判断出来了），此时需要看轨道图形来判断
作者Author: 我本是个娃娃 时间: 2015-6-28 11:20
觉得不用来路不明的关键词应该是遵循奥卡姆剃刀原理吧
作者Author: limaolin0 时间: 2015-6-29 23:26
又收获一点点了
作者Author: lastzealot 时间: 2015-6-30 21:30
写得好 受教了
作者Author: beefly 时间: 2015-7-2 02:05
本帖最后由 beefly 于 2015-7-2 02:14 编辑
用nosymm做频率计算还会导致一个危险，就是热化学量（熵、自由能）是错误的，因为它们依赖转动熵，而转动熵与点群有关。最极端的情况下，用和不用对称性，自由能相差好几个kcal/mol。例如水分子，标准状态下的自由能相差0.4 kcal/mol。C60：2.5 kcal/mol。CO2：2.4 kcal/mol。用了比实际低的对称性，也有同样的问题。如果是很高的温度，误差就更大了。
但是Gaussian频率部分的对称性程序比较独立，并不完全受nosymm的影响。在热化学开始的部分，找到一行
Rotational symmetry number ...
从后面的数值可以判断热化学计算是否用了正确的对称性。
对于线形分子，在优化的时候如果nosymm导致线形结构被破坏，那么频率和热化学计算的时候除了以上错误以外，还会导致丢掉一个振动模式和相应的各种热化学贡献。
作者Author: 神龍 时间: 2015-7-26 21:15
感谢sob老师！
作者Author: GoldenBaby 时间: 2019-10-1 14:37
讲一个比较特殊的例子，有些及其特殊的情况，利用对称性容易导致scf不收敛，这种情况下使用nosymm可以得到解决（PS：因为最近在做一个类似结构的优化的时候发现在该计算级别下，在sob老师在博文里提到的帮助scf收敛的方法，无论用哪种方法都无法收敛，但是事实上在更低的计算级别下可以收敛的，遂尝试使用nosymm关键词，发现可以收敛并最终完成优化，因此讲这个例子分享给大家，同时如果各位有什么更好的意见或建议欢迎批评指正。）
作者Author: zjxitcc 时间: 2019-10-2 14:32
本帖最后由 zjxitcc 于 2019-10-2 14:40 编辑
看了你的文件，有些建议，不过与本帖主题nosymm无关。你这个计算虽然收敛了，但是RHF能量并非最低，且所有的RHF波函数都有RHF -> UHF不稳定性。铁钴镍这些过渡金属基本上都要用UDFT来描述的。特别是你算的这个还是低自旋（对[Fe(CN)6]4-而言，低自旋能量是比较高的，具有明显多组态特征）。在你这个坐标、方法和基组下有几个SCF解：
method Elec E <S**2> stability
RPBE0 -680.082798 0 RHF->UHF instability
RPBE0 -680.091368 0 RHF->UHF instability
UPBE0 -680.091884 0 internal instability
UPBE0 -680.091893 0.8394 internal instability
UPBE0 -680.146410 1.0157 stable
UPBE0 -680.147835 1.0284 stable
如果是研究[Fe(CN)6]4-不同自旋态下的能量差，高自旋你用UDFT算，低自旋用RDFT算显然不合适，这比较不公平；一旦低自旋也用UDFT，那么UDFT同样有多个SCF解，通常你只能取能量最低的解。
比较严谨的方式是在基态（一般是高自旋）下就用CASSCF算，然后其他自旋用CASSCF的不同自旋多重度的root来考虑，这样每个自旋多重度都可以考虑到多组态特征。如果不想做/不会做/算不动CASSCF的话，就用UDFT也凑合，就是要注意SCF解有好几个。
PS：这个例子从一定程度上说明了scf=qc虽然收敛性好，但是由于其二阶性质容易收敛到最近的local minimum上，并不一定是全局极小。
附上我的计算文件 链接: https://pan.baidu.com/s/1jIYEOLaUHX4U3Yfj6Lql6A 提取码: 746q
作者Author: GoldenBaby 时间: 2019-10-3 00:35
感谢大神分享，学习了
作者Author: sobereva 时间: 2019-10-3 04:18
和nosymm没有直接关系。单纯地想让限制性闭壳层下的PBE0/SDD_6-311+G*的SCF收敛，g16依次运行以下两个文件就能实现，没利用我强烈不主张优先考虑的scf=qc，也没有利用nosymm
(, 下载次数 Times of downloads: 65)
(, 下载次数 Times of downloads: 61)
但如第二个任务对波函数稳定性的测试，当前存在波函数不稳定性
如果把第二个任务加上stable=opt，会自动给出稳定的自旋极化单重态波函数
(, 下载次数 Times of downloads: 29)
有的人表面上看加了nosymm就收敛了，就以为nosymm对帮助SCF收敛有益、应当作为解决SCF不收敛应当要考虑的常规做法中，这是严重的误解。
作者Author: GoldenBaby 时间: 2019-10-3 14:51
感谢sob老师的回复，学习了
作者Author: ionexchangeC 时间: 2022-4-29 00:40
这里想补充一点：对于G16，如果想在结构优化过程中保持对称性不降低，还有一种方法是利用广义内坐标GIC。大致的方法是，分析Gaussian结构优化开始时建立的冗余内坐标，指出在欲保持的点群下哪些变量是相等的，哪些是相加等于180度或者360度的，然后设置一些冻结的广义内坐标，保证这些数学关系的成立。
如：令B(1,2)和B(1,4)两个键长变量相等，可设置F1(freeze)=B(1,2)-B(1,4)，并使他们在输入文件中就相等，优化时因为差值固定所以保持相等。
令A(1,2,3)、A(3,2,4)和A(4,2,1)相加等于360度，可设置F2(freeze)=A(1,2,3)+A(3,2,4)+A(4,2,1)，并在输入文件中保证三键角相加等于360度。
作者Author: sobereva 时间: 2022-4-29 07:57
用GIC实现往往遇到的一个实际问题是优化中途容易莫名其妙出错，和约束具体定义方式很有关系
作者Author: gauss98 时间: 2022-12-21 12:32
“用nosymm做频率计算还会导致一个危险，就是热化学量（熵、自由能）是错误的，”
这个对较复杂在C1对称性结构也是这样的吗？
我最近发现用不用 nosymm 对自由能计算有差异（C1对称性结构）
这么说来应该是不用nosymm的结果更可信？
作者Author: JCenter 时间: 2024-11-5 11:25
老师，我在计算电子密度差时，使用了b3lyp-d3(bj)/def2TZVP结合iefpcm模型代表水环境的计算级别得到波函数信息，而且使用了nosymm。想问下，nosymm会影响体系的电子密度分布吗，想直接使用已有的波函数文件基于Multiwfn接着计算ADCH原子电荷，不知是否可以呢。
作者Author: sobereva 时间: 2024-11-8 21:48
要注意看热力学计算阶段对点群的判断，Gaussian直接输出了。不管用不用nosymm，那个阶段都会判断点群，从而试图恰当确定转动对称数，以计算转动对热力学量的贡献。
(, 下载次数 Times of downloads: 0)
作者Author: sobereva 时间: 2024-11-8 21:50
没直接影响
可以
作者Author: JCenter 时间: 2024-11-8 22:41
感谢老师
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |