一个用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
    

复制代码

wzkchem5

近两年ORCA可能会实现DLPNO-CCSD(T)解析梯度，我们这边有一个博士生正在写。在此之前会实现canonical CCSD(T)解析梯度。
此外有兴趣的话可以看看NEWUOA和ORCA内置的基于数值梯度的结构优化算法哪个快。对于18碳环这样的高对称性分子肯定前者快，但是可以测一下非对称性分子，假如也是NEWUOA快的话，可以考虑写一个基于NEWUOA 的通用结构优化脚本，上传到这里https://gitlab.gwdg.de/orca-helpers/orca-helpers/-/tree/master/（这里有很多普通orca用户写的脚本），这样其他人也可以用

sobereva

顺带一提，一些范围分离双杂化泛函可以正确优化C18，如SOS1-QIDH和RSX-QIDH

hebrewsnabla

对于DLPNO系列方法，包括DLPNO-MP2和DLPNO-CC，数值梯度可能有比较大的问题，参见 https://orcaforum.kofo.mpg.de/vi ... 7610&p=32831#p32798 ，原因是数值扰动时如果轨道domain改变，势能面就不够光滑。而解析梯度（如果存在的话）不会有这个问题。

虽然NEWUOA不是基于数值梯度的，但是搞不好也有这个问题。

mfdsrax2

wzkchem5 发表于 2022-7-14 01:34
近两年ORCA可能会实现DLPNO-CCSD(T)解析梯度，我们这边有一个博士生正在写。在此之前会实现canonical CCSD( ...

这个要写两年吗

zjxitcc

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

两年能写出来就已经很好了，这巨复杂。比如MOLCAS的CASPT2解析导数做了N年也没做出来

mfdsrax2

zjxitcc 发表于 2022-7-14 10:52
两年能写出来就已经很好了，这巨复杂。比如MOLCAS的CASPT2解析导数做了N年也没做出来

DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂

zjxitcc

mfdsrax2 发表于 2022-7-14 10:54
DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂

即使是正则CCSD(T)解析导数也很复杂。而局域相关方法的解析导数复杂得多，即使写出来了，还要考虑解析导数是不是也实现了线性标度的问题。

mfdsrax2

zjxitcc 发表于 2022-7-14 11:00
即使是正则CCSD(T)解析导数也很复杂。而局域相关方法的解析导数复杂得多，即使写出来了，还要考虑解析导 ...

那为什么不先尝试CCSD(T)的解析导数，而是花两年做更难的DLPNO-CCSD(T)导数。CCSD(T)的解析导数对于计算化学的意义也很大

zjxitcc

mfdsrax2 发表于 2022-7-14 11:07
那为什么不先尝试CCSD(T)的解析导数，而是花两年做更难的DLPNO-CCSD(T)导数。CCSD(T)的解析导数对于计算 ...

2L里有说“在此之前会实现canonical CCSD(T)解析梯度”，这就是正则CCSD(T)，或者说传统CCSD(T)。传统CCSD(T)解析导数已经有好几个量化程序做出来了，ORCA实现它只能是增添一个小功能（因为具有可替代性，且传统CCSD(T)标度是N^7，原理上就没法算大体系），作为自己写代码时的验证。重要的是DLPNO-CCSD(T)解析导数。

wzkchem5

mfdsrax2 发表于 2022-7-14 03:54
DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂

基本上任何近似方法的代码都比相应的非近似方法的代码复杂。至少当你只需要做对而不需要做快的时候是如此的。

Mikasa

zjxitcc 发表于 2022-7-14 10:52
两年能写出来就已经很好了，这巨复杂。比如MOLCAS的CASPT2解析导数做了N年也没做出来

现在的OpenMOLCAS有CASPT2的解析梯度了呀

biogon

Mikasa 发表于 2022-7-15 09:10
现在的OpenMOLCAS有CASPT2的解析梯度了呀

我看22jul08的手册还说没有啊

biogon

wzkchem5 发表于 2022-7-14 17:00
基本上任何近似方法的代码都比相应的非近似方法的代码复杂。至少当你只需要做对而不需要做快的时候是如此 ...

orca的NEVPT2梯度现在有没有人在做，以前论坛上说打算做了

Mikasa

biogon 发表于 2022-7-15 09:39
我看22jul08的手册还说没有啊

可能只是手册没更新，我们已经用上了~