一个用DLPNO-CCSD(T)优化18碳环的结构的思路

ionexchangeC

由于DLPNO-CCSD(T)没有解析梯度，常规的优化方式没法用，ORCA又不能利用对称性在内坐标下优化。于是我自己写了一个调用ORCA的小程序实现D9h结构的18碳环的结构优化。由于DLPNO-CCSD(T)结合像样的基组在我的电脑上算个单点都很花时间，我换成了PWPB95-D4/def2-TZVP这个级别（也没有解析梯度）。优化算法是NEWUOA。
18C.f90 (51.76 KB, 下载次数 Times of downloads: 9)

2022-7-13 23:07 上传 Uploaded

点击下载
Click to download

18C.exe (434.87 KB, 下载次数 Times of downloads: 0)

2022-7-13 23:07 上传 Uploaded

点击下载
Click to download

当然最后优化的结果是碳碳键键长几乎相等，是不合理的。另外，由于18碳环的HF波函数定性错误（RHF有波函数不稳定问题），我觉得可以用DFT波函数作为参考态，应当是更合理的。ORCA直接就支持（写M062X DLPNO-CCSD(T)之类）。
程序的核心代码如下：

1. PROGRAM test_newuoa
   
2. use newuoa_module
   
3. implicit real*8 (a-h,o-z)
   
4. real*8,allocatable :: X(:)
   
5. external CALFUN
   
6. 
   
7. nvar=2
   
8. allocate(X(nvar))
   
9. X(:)=0
   
10. X(1)=1.3436D0 !长碳碳键键长初猜
    
11. X(2)=1.2242D0 !短碳碳键键长初猜
    
12. call newuoa_min(CALFUN, X, RHOBEG=0.1D-1, RHOEND=1D-4, IPRINT=2, MAXFUN=50000) !键长变化一般不会超过0.1A，所以RHOBEG设置了0.1D-1
    
13. END PROGRAM
    
14.         
    
15. subroutine CALFUN(X,F)
    
16. real*8 :: X(:),F
    
17. real*8 ene
    
18. open(168,FILE="./4.inp")
    
19. write(168,*) "%pal nprocs 4 end"
    
20. write(168,*) ""
    
21. write(168,*) "! PWPB95 D4 def2-TZVP RIJK def2/JK def2-TZVP/C TightSCF" !计算级别行
    
22. write(168,*) "%maxcore 1000"
    
23. write(168,*) "* gzmt 0 1"
    
24. write(168,*) " C              "
    
25. write(168,*) " C                  1            ",X(1)
    
26. write(168,*) " C                  2            ",X(2),"    1            160.00000000"
    
27. write(168,*) " C                  3            ",X(1),"    2            160.00000000    1            0.00000000    "
    
28. write(168,*) " C                  4            ",X(2),"    3            160.00000000    2            0.00000000    "
    
29. write(168,*) " C                  5            ",X(1),"    4            160.00000000    3            0.00000000    "
    
30. write(168,*) " C                  6            ",X(2),"    5            160.00000000    4            0.00000000    "
    
31. write(168,*) " C                  7            ",X(1),"    6            160.00000000    5            0.00000000    "
    
32. write(168,*) " C                  8            ",X(2),"    7            160.00000000    6            0.00000000    "
    
33. write(168,*) " C                  9            ",X(1),"    8            160.00000000    7            0.00000000    "
    
34. write(168,*) " C                 10           ",X(2),"    9            160.00000000    8            0.00000000    "
    
35. write(168,*) " C                 11           ",X(1),"   10           160.00000000    9            0.00000000    "
    
36. write(168,*) " C                 12           ",X(2),"   11           160.00000000   10           0.00000000    "
    
37. write(168,*) " C                 13           ",X(1),"   12           160.00000000   11           0.00000000    "
    
38. write(168,*) " C                 14           ",X(2),"   13           160.00000000   12           0.00000000    "
    
39. write(168,*) " C                 15           ",X(1),"   14           160.00000000   13           0.00000000    "
    
40. write(168,*) " C                 16           ",X(2),"   15           160.00000000   14           0.00000000    "
    
41. write(168,*) " C                 17           ",X(1),"   16           160.00000000   15           0.00000000    " !内坐标
    
42. write(168,*) "*"
    
43. close(168)
    
44. call system('E:\Baidunetdiskdownload\orca503\orca 4.inp > 4.out') !运行ORCA，路径写ORCA路径（我是windows下运行的）
    
45. call system("grep 'FINAL SINGLE' 4.out | awk '{print $5}' > 4.txt") !提取数据
    
46. open(10,file="4.txt")
    
47. read(10,*) ene
    
48. print "(3D20.12)",ene,X(1),X(2) !输出数据
    
49. close(10)
    
50. call system('rm 4.txt')
    
51. F=ene
    
52. end subroutine

复制代码

最后程序的输出是

1. -0.685090685940D+03  0.134360000000D+01  0.122420000000D+01
   
2. -0.685088824064D+03  0.135360000000D+01  0.122420000000D+01
   
3. -0.685091935687D+03  0.134360000000D+01  0.123420000000D+01
   
4. -0.685091376254D+03  0.133360000000D+01  0.122420000000D+01
   
5. -0.685086911391D+03  0.134360000000D+01  0.121420000000D+01
   
6. -0.685093443649D+03  0.133360049157D+01  0.123410084787D+01
   
7. -0.685094798943D+03  0.132032800904D+01  0.124176231458D+01
   
8. -0.685094995375D+03  0.131259737079D+01  0.124433184094D+01
   
9. -0.685095034243D+03  0.131324281292D+01  0.125431098942D+01
   
10. -0.685095381505D+03  0.130901412440D+01  0.125156104932D+01
    
11. -0.685095043825D+03  0.130676166617D+01  0.126130406897D+01
    
12. -0.685095617327D+03  0.130113934072D+01  0.125772447223D+01
    
13. -0.685096084259D+03  0.128547177783D+01  0.126955992432D+01
    
14. -0.685094926437D+03  0.125415034613D+01  0.129324894559D+01
    
15. -0.685094328761D+03  0.126621797661D+01  0.127341240296D+01
    
16. -0.685095912811D+03  0.128405212975D+01  0.127945864137D+01
    
17. -0.685096232455D+03  0.127922217111D+01  0.127736648668D+01
    
18. -0.685096136254D+03  0.127312956531D+01  0.128529618747D+01
    
19. 
    
20.     New RHO = 1.0000D-03     Current number of function evaluations =    18
    
21.     Least value of F = -6.850962324546280D+02
    
22.     The corresponding X array is:
    
23.      1.279222D+00   1.277366D+00
    
24. -0.685096233410D+03  0.127758652046D+01  0.127882988316D+01
    
25. -0.685096234289D+03  0.127809725865D+01  0.127837670778D+01
    
26. -0.685096220707D+03  0.127736049474D+01  0.127770055769D+01
    
27. -0.685096223664D+03  0.127897725857D+01  0.127885168162D+01
    
28. 
    
29.     New RHO = 1.0000D-04     Current number of function evaluations =    22
    
30.     Least value of F = -6.850962342889070D+02
    
31.     The corresponding X array is:
    
32.      1.278097D+00   1.278377D+00
    
33. -0.685096234359D+03  0.127830643932D+01  0.127823790676D+01
    
34. -0.685096234358D+03  0.127837923057D+01  0.127816933975D+01
    
35. -0.685096234369D+03  0.127823904017D+01  0.127816403285D+01
    
36. -0.685096233929D+03  0.127816395443D+01  0.127809798641D+01
    
37. -0.685096234265D+03  0.127830070069D+01  0.127817211338D+01
    
38. -0.685096234304D+03  0.127816288187D+01  0.127822883960D+01
    
39. -0.685096234345D+03  0.127830944408D+01  0.127816255353D+01
    
40. 
    
41.     At the return from NEWUOA     Total times of function evaluations =    29
    
42.     Least value of F = -6.850962343688360D+02
    
43.     The corresponding X array is:
    
44.      1.278239D+00   1.278164D+00
    

复制代码

Eclair

mfdsrax2 发表于 2022-7-14 10:34
这个要写两年吗

你觉得慢可以自己来写，ccsd（T）的解析梯度很早以前就推出来了，

ionexchangeC

chrinide 发表于 2023-12-14 23:20
这是算了多长时间了？ 不过结果可以用来引用的

我是为了测试超算做的这个任务，而我现在甚至忘记了测试账号的密码

chrinide

ionexchangeC 发表于 2023-12-12 19:35
2023.12.12更新：
用RHF为参考态，DLPNO-CCSD(T)/cc-pVTZ级别优化18碳环，结果是长键1.359埃，短键1.232埃

这是算了多长时间了？ 不过结果可以用来引用的

ionexchangeC

2023.12.12更新：
用RHF为参考态，DLPNO-CCSD(T)/cc-pVTZ级别优化18碳环，结果是长键1.359埃，短键1.232埃

Kalinite

Mikasa 发表于 2022-7-15 09:43
可能只是手册没更新，我们已经用上了~

您好，请问应该通过什么方式实现呢，我在输入里直接调ALASKA，输出文件中显示仍然是numerical gradient。

1. Numerical gradient, root     1
   
2.   ---------------------------------------------
   
3.                  X           Y           Z
   
4.   ---------------------------------------------
   
5.   O1          -0.046735   -0.066128   -0.000000
   
6.   H1           0.086732   -0.011718   -0.000000
   
7.   H2          -0.039997    0.077846    0.000000
   
8.   ---------------------------------------------
   
9. --- Stop Module: numerical_gradient at Mon Dec 19 10:04:31 2022 /rc=_RC_ALL_IS_WELL_ ---
   
10. *** files: xmldump
    
11.     saved to directory /export/home/chli/xjw/test/molcas/caspt2grad
    
12. --- Module numerical_gradient spent 17 seconds ---
    

复制代码

wjc404

ionexchangeC 发表于 2022-11-8 11:39
已上传CASSCF波函数，欢迎讨论

谢谢。我用ORCA跑cas(8,8)出来的波函数和这个差不多，能够收敛，但是后续的MCRPA计算发现casscf波函数不稳定（以及nevpt2计算中警告Koopmans矩阵存在负的本征值（结果未上传））。
计算结果: https://pan.baidu.com/s/1mL_uUeRGpgR_oknlpXv9rg 提取码: ari5

ionexchangeC

wjc404 发表于 2022-11-7 23:47
看了一下（R1=1.22 Å，R2=1.34 Å时） RHF/def2TZVP的轨道hessian的（共3个）本征值小于零对应 ...

已上传CASSCF波函数，欢迎讨论

ionexchangeC

cyclo[18]carbon.molden.input (7.48 MB, 下载次数 Times of downloads: 3)

2022-11-8 11:38 上传 Uploaded

点击下载
Click to download

wjc404

ionexchangeC 发表于 2022-8-15 11:48
又用CASSCF和NEVPT2优化了18碳环，取的活性空间是（8，8）。
CASSCF：长键1.356埃，短键1.199埃
NEVPT2： ...

看了一下（R1=1.22 Å，R2=1.34 Å时） RHF/def2TZVP的轨道hessian的（共3个）本征值小于零对应的向量中大于0.1的元素，发现牵涉最多的轨道为最高的8个占有轨道（均为π）和最低的8个未占有π轨道。它们在每个向量中都有较高贡献，且互相之间没有绝对的主次。猜测在轨道优化时仅考虑（8，8）的活性空间可能偏小了。但再往大就比较昂贵了。

ionexchangeC

wjc404 发表于 2022-10-25 11:44
是否可以用CCSD(T)-F12(RI)/cc-pVDZ-F12试一下优化(看情况启用DLPNO；这个级别的精度相当于不带F12时的TZ基 ...

按理来说，你这个计算级别就算加了DLPNO也比普通地用DLPNO-CCSD(T)/cc-pVTZ贵吧

wjc404

是否可以用CCSD(T)-F12(RI)/cc-pVDZ-F12试一下优化(看情况启用DLPNO；这个级别的精度相当于不带F12时的TZ基组结果)？

另外想请教一下，现在有程序支持“DFT的基态上做CC方法”的解析梯度吗？

ionexchangeC

又用CASSCF和NEVPT2优化了18碳环，取的活性空间是（8，8）。
CASSCF：长键1.356埃，短键1.199埃
NEVPT2：长键等于短键，1.289埃
CASSCF.txt (2.34 KB, 下载次数 Times of downloads: 4)

2022-8-15 11:48 上传 Uploaded

点击下载
Click to download

FIC-NEVPT2.txt (2.58 KB, 下载次数 Times of downloads: 7)

2022-8-15 11:48 上传 Uploaded

点击下载
Click to download

sobereva

zjxitcc 发表于 2022-8-10 23:09
为啥不用基于自旋极化UHF的DLPNO-UCCSD(T)。。

DLPNO-UCCSD(T)比DLPNO-CCSD(T)昂贵甚巨
而且UHF错误地高估了这个体系的对称破缺特征，用U还不如用R
另外U和R我实际用CCSD优化的结果基本没差别

zjxitcc

ionexchangeC 发表于 2022-8-10 22:39
用DLPNO-CCSD(T)-M06-2X/cc-pVTZ（M06-2X为参考态的DLPNO-CCSD(T)）按照这个方法优化，最后的结果是长碳碳 ...

为啥不用基于自旋极化UHF的DLPNO-UCCSD(T)。。