请问Gaussian16可以计算基于KS轨道的CCSD能量吗

wjc404

发现ORCA, Q-CHEM和MRCC是可以计算基于KS轨道的CCSD能量，而且相互符合得很好。
比如在计算NO2(R=1.19, A=134.0)的基态时：
ORCA写如下关键词

1. UKS B3LYP/G CCSD(T) def2-TZVP ExtremeSCF UseSym DefGrid3 NoRI NoFrozenCore

复制代码

Q-CHEM串联两个任务分别写如下关键词

1. $rem
   
2.    JOBTYPE sp
   
3.    THRESH 14
   
4.    METHOD b3lyp
   
5.    XC_GRID 000250000974
   
6.    SCF_CONVERGENCE 8
   
7.    STABILITY_ANALYSIS TRUE
   
8.    BASIS def2-TZVP
   
9. $end

复制代码

1. $rem
   
2.    JOBTYPE sp
   
3.    SCF_GUESS read
   
4.    METHOD ccsd
   
5.    CC_CONVERGENCE 7
   
6.    MAX_SCF_CYCLES 0
   
7.    CC_CANONIZE TRUE
   
8.    N_FROZEN_CORE 0
   
9.    BASIS def2-TZVP
   
10. $end

复制代码

MRCC写如下关键词

1. basis=def2-TZVP
   
2. calc=CCSD(T)
   
3. dft=B3LYP3
   
4. agrid=LD0110-LD0974
   
5. grtol=12
   
6. scftol=9
   
7. ccprog=mrcc
   
8. core=corr
   
9. cctol=8
   
10. dfbasis_scf=none
    
11. dfbasis_cor=none

复制代码

三者得出的能量，B3LYP的均为-205.167120 Hartree，CCSD的均为-204.825489 Hartree。（CCSD(T)的小数点后第三位不同，这里暂不讨论）
然后尝试用Gaussian16做，先计算KS轨道（B3LYP能量与其他程序一致），再在下一个任务读取轨道算后HF相关。在第二个任务中，用nonstd的方式写以下关键词

1. 1/29=7,38=1/1;
   
2. 2/12=2,40=1/2;
   
3. 3/5=7,6=2,11=9,25=1,30=1,116=-2/1,2,3;
   
4. 4/5=1,9=1000,117=-1/1;
   
5. 8/6=4,9=120000,10=90/1,4;
   
6. 9/5=7,9=9,14=2/13;
   
7. 99/5=1,9=1/99;

复制代码

发现结果中reference determinant的能量是和ORCA给出的一样的，但CCSD能量-204.8289 Hartree就有些问题（而且最大振幅是0.102，不同于ORCA和MRCC的0.088）。
然后尝试做基于B3LYP优化的轨道的CCD计算，Gaussian16的第二个任务输入中nonstd的第六行改为"9/5=6,9=9/13;"，ORCA的任务将关键词CCSD(T)改为CCD。发现两者算出的CCD能量均为-204.817902 Hartree。
另外发现如果以上Gaussian16的计算中第一步是UHF而不是UKS，则第二步算出的CCSD能量与单步直接算UCCSD是一样的。


从以上结果看，应该是Gaussian16计算CCSD时non-HF开关没打开（或者没有实现相应的代码），然后从IOp文档(l801.hlp, l913.hlp)中也没看到应该从哪里设置（试了比较相近的iop(8/46=1)发现没用）。想请教一下对于Gaussian有没有正确的方式能做KS轨道的CCSD能量？

niobium

你用ORCA和MRCC都实现了，还折腾Gaussian干啥呢

wjc404

niobium 发表于 2023-1-7 20:20
你用ORCA和MRCC都实现了，还折腾Gaussian干啥呢

谢谢关注。主要是这两个程序算稍大的体系有点效率问题。ORCA基于MPI并行时不同进程的内存共享不多，很多积分内存放不下要放到硬盘上，耦合簇迭代时硬盘IO稍夸张，然后DLPNO不适合自旋极化单重态。MRCC的代码不是为CCSD专门优化的（有一个专门做CCSD的程序，但写了只支持闭壳层以及不支持对称性），时间上会慢一个数量级。Gaussian16的后HF积分变换及CCSD代码的效率还是相对比较高的，内存和硬盘使用上也比较稳健。

zjxitcc

用DFT轨道做CC计算是邪典，强烈建议所有量化工作者放弃此做法，以免走入歧途（除非是想证明它结果不好 而去学习使用它）。体系如果有多参考态特征，应该采用多组态/多参考态方法。即使是单参考方法，想要算更大的体系，还有DLPNO-CCSD(T)；想要更高精度，还有CCSD(T)-F12。

看《详细对比基于HF和KS-DFT轨道做耦合簇的一篇文章（JCC,2022）》http://bbs.keinsci.com/thread-32791-1-1.html，其中提到“基于KS-DFT轨道做耦合簇并不会带来任何好处，有时结果更好也是由于误差抵消带来的。仅当HF计算遇到SCF或耦合簇方程迭代不收敛的时候，才值得尝试基于KS-DFT轨道做耦合簇”。

wjc404

zjxitcc 发表于 2023-1-8 16:04
用DFT轨道做CC计算是邪典，强烈建议所有量化工作者放弃此做法，以免走入歧途。体系如果有多参考态特征，应 ...

谢谢回复。我只是发现Gaussian的CCSD没有显式支持 non-HF reference 的情况（此时方程的某些和单激发相关的展开项不为零，可能程序默认它们是零了），在此贴中以KS轨道来举例。KS轨道在理论上作为CC的reference确实比较奇怪，两者八竿子打不到一块。
有时候HF作为CC的reference会有问题，比如一些 artificial symmetry breaking 的体系。Gaussian有Brueckner-CCD可以处理这些体系的情况，但有时会出现收敛巨慢的问题（比如每次轨道优化的结果是residual乘了0.8，则要轨道优化10次以上才能让残差降一个数量级！然后每次轨道优化都涉及积分变换的重算及若干次双激发幅度的迭代优化），导致时间的消耗比同基组的CCSD高两个数量级。。

喵星大佬

NWChem也可以，效率其实也还行，但这种歪门邪道嘛

zjxitcc

用MOKIT自动做NEVPT2/def2TZVP计算，16核耗时仅18s，我相信比CC计算快不少。计算过程中自动确定活性空间是(3,3)，CASSCF自然轨道及其占据数如下


MOKIT输入文件NO2.gjf内容如下

1. %mem=16GB
   
2. %nprocshared=16
   
3. #p NEVPT2/def2TZVP
   
4. 
   
5. mokit{}
   
6. 
   
7. 0 2
   
8. N    -1.14809785    0.42119565    0.0
   
9. O     0.00655407    0.13330859    0.0
   
10. O    -1.74309785    1.45176588    0.0

复制代码

E(UHF) =      -204.12100456 a.u., <S**2>=  0.770
E(CASSCF) =      -204.14340332 a.u.
E(SC-NEVPT2) =      -204.83635317 a.u. （无冻芯）
可供参考和对比。

wzkchem5

zjxitcc 发表于 2023-1-8 09:04
用DFT轨道做CC计算是邪典，强烈建议所有量化工作者放弃此做法，以免走入歧途。体系如果有多参考态特征，应 ...

我同意当HF determinant定性正确时，KS reference统计上不会明显比HF reference好。但是有的体系是HF reference定性错误，而KS reference定性正确，此时KS reference可能还是有帮助的，问题可能仅在于（1）这类体系比较少；（2）这类体系一般有显著的多参考性质（但改用多参考态方法未必现实，第一可能算不动，第二动态相关不可避免地会处理得更粗糙）

zjxitcc

wzkchem5 发表于 2023-1-9 18:09
我同意当HF determinant定性正确时，KS reference统计上不会明显比HF reference好。但是有的体系是HF ref ...

我想知道“HF reference定性错误，而KS reference定性正确”有啥例子我认为几乎不可能存在，因为复杂体系UHF有多重解，可能得出这个结论的人用了明显错误的UHF/UDFT解（例如波函数不稳定）而不自知。不过这也只是我的猜测，期待有一天有这种例子可以拿来具体算一算。

beefly

zjxitcc 发表于 2023-1-9 21:03
我想知道“HF reference定性错误，而KS reference定性正确”有啥例子我认为几乎不可能存在，因为复 ...

很多高对称的开壳层体系做HF计算都存在artifactual broken symmetry问题。以三原子ABA类型的开壳层体系为例（为了简单假设是两个AB键长相等的直线构型，但这些条件不是必须的），基态波函应当是AB-A和A-BA两种组态的线性组合，各占50%。在进行HF计算时，当把D2h（或Dooh）对称性降到C2v（或Coov）会发现，HF只能收敛到其中一个组态，C2v对称性的HF能量比D2h对称性的HF能量更低。前者是HF缺少电子关联导致的非物理解（例如，会出现超大的电子偶极矩和振动频率），因此后续的MP2、CC也没有意义。虽然MCSCF原则上可以解决此问题，但是必须充分考虑静态关联，否则仍可能得到非物理解。而考虑了电子关联的DFT则很少出现这种问题，个别体系用杂化因子比较高的杂化泛函可以避免。

代表体系：FHF，BNB，ZnCl2+，CO2和CuCl2的某些激发态。molcas的论坛上曾经讨论过的O4也属于这种情况，不过现在看不到了

scf

有些DMRG计算貌似用的BP86的轨道， 比如https://www.nature.com/articles/nchem.2041#Sec7的SI，"To generate the active space for the DMRG calculations, we performed an unrestricted DFT BP86/SVP calculation for the high spin (Sz =5) state. The alpha occupied and unoccupied or-
bitals were then separately localized (“split-localized”) [4] using the Pipek-Mezey algorithm [5]. "不知道为什么没用HF轨道

wjc404

wjc404 发表于 2023-1-9 13:53
谢谢回复。我只是发现Gaussian的CCSD没有显式支持 non-HF reference 的情况（此时方程的某些和单激发相关 ...

发现Gaussian16的ILSW（一个长度为200且每个元素8字节的控制向量）的第22个元素像是一个与non-HF相关的开关，但似乎光把它开启还不够，还要准备752文件（以MO为基的Fock矩阵，alpha和beta各有一份，各自的只保留占有轨道（不包括冻芯）和虚轨道之间的（对称性允许的）耦合系数），以及可能还需要760文件（semi-canonicalize后的MO，alpha和beta各一份）。这两个文件应该是l801.exe识别到波函数为ROHF后去自动产生的，其他方式看起来行不通。程序写死了，就没有办法。752文件中每个占有轨道的关联矩阵元都挨在一起。占有轨道之间的排序优先级为（前高后低）：alpha先beta后，不可约表示小的先，能量低的先。

zjxitcc

wjc404 发表于 2023-1-11 23:50
发现Gaussian16的ILSW（一个长度为200且每个元素8字节的控制向量）的第22个元素像是一个与non-HF相关的开 ...

我也构思过大概做法，没你这么详细，想到最后觉得不如自己写一个CCSD程序