计算化学公社

标题: 一个用DLPNO-CCSD(T)优化18碳环的结构的思路 [打印本页]
作者
Author: ionexchangeC    时间: 2022-7-13 23:29
标题: 一个用DLPNO-CCSD(T)优化18碳环的结构的思路
由于DLPNO-CCSD(T)没有解析梯度，常规的优化方式没法用，ORCA又不能利用对称性在内坐标下优化。于是我自己写了一个调用ORCA的小程序实现D9h结构的18碳环的结构优化。由于DLPNO-CCSD(T)结合像样的基组在我的电脑上算个单点都很花时间，我换成了PWPB95-D4/def2-TZVP这个级别（也没有解析梯度）。优化算法是NEWUOA。
(, 下载次数 Times of downloads: 9) (, 下载次数 Times of downloads: 0)
当然最后优化的结果是碳碳键键长几乎相等，是不合理的。另外，由于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

  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

  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

  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

  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
复制代码
作者
Author: wzkchem5    时间: 2022-7-14 01:34
近两年ORCA可能会实现DLPNO-CCSD(T)解析梯度，我们这边有一个博士生正在写。在此之前会实现canonical CCSD(T)解析梯度。
此外有兴趣的话可以看看NEWUOA和ORCA内置的基于数值梯度的结构优化算法哪个快。对于18碳环这样的高对称性分子肯定前者快，但是可以测一下非对称性分子，假如也是NEWUOA快的话，可以考虑写一个基于NEWUOA 的通用结构优化脚本，上传到这里https://gitlab.gwdg.de/orca-helpers/orca-helpers/-/tree/master/（这里有很多普通orca用户写的脚本），这样其他人也可以用
作者
Author: sobereva    时间: 2022-7-14 06:29
顺带一提，一些范围分离双杂化泛函可以正确优化C18，如SOS1-QIDH和RSX-QIDH
作者
Author: hebrewsnabla    时间: 2022-7-14 09:56
本帖最后由 hebrewsnabla 于 2022-7-14 10:06 编辑

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

虽然NEWUOA不是基于数值梯度的，但是搞不好也有这个问题。
作者
Author: mfdsrax2    时间: 2022-7-14 10:34
wzkchem5 发表于 2022-7-14 01:34
近两年ORCA可能会实现DLPNO-CCSD(T)解析梯度，我们这边有一个博士生正在写。在此之前会实现canonical CCSD( ...

这个要写两年吗
作者
Author: zjxitcc    时间: 2022-7-14 10:52
mfdsrax2 发表于 2022-7-14 10:34
这个要写两年吗

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

DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂
作者
Author: zjxitcc    时间: 2022-7-14 11:00
mfdsrax2 发表于 2022-7-14 10:54
DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂

即使是正则CCSD(T)解析导数也很复杂。而局域相关方法的解析导数复杂得多，即使写出来了，还要考虑解析导数是不是也实现了线性标度的问题。
作者
Author: mfdsrax2    时间: 2022-7-14 11:07
zjxitcc 发表于 2022-7-14 11:00
即使是正则CCSD(T)解析导数也很复杂。而局域相关方法的解析导数复杂得多，即使写出来了，还要考虑解析导 ...

那为什么不先尝试CCSD(T)的解析导数，而是花两年做更难的DLPNO-CCSD(T)导数。CCSD(T)的解析导数对于计算化学的意义也很大
作者
Author: zjxitcc    时间: 2022-7-14 11:11
本帖最后由 zjxitcc 于 2022-7-14 11:13 编辑
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)解析导数。
作者
Author: wzkchem5    时间: 2022-7-14 17:00
mfdsrax2 发表于 2022-7-14 03:54
DLPNO-CCSD(T)相比CCSD(T)的解析导数是更容易还是复杂

基本上任何近似方法的代码都比相应的非近似方法的代码复杂。至少当你只需要做对而不需要做快的时候是如此的。
作者
Author: Mikasa    时间: 2022-7-15 09:10
zjxitcc 发表于 2022-7-14 10:52
两年能写出来就已经很好了，这巨复杂。比如MOLCAS的CASPT2解析导数做了N年也没做出来

现在的OpenMOLCAS有CASPT2的解析梯度了呀
作者
Author: biogon    时间: 2022-7-15 09:39
Mikasa 发表于 2022-7-15 09:10
现在的OpenMOLCAS有CASPT2的解析梯度了呀

我看22jul08的手册还说没有啊
作者
Author: biogon    时间: 2022-7-15 09:40
wzkchem5 发表于 2022-7-14 17:00
基本上任何近似方法的代码都比相应的非近似方法的代码复杂。至少当你只需要做对而不需要做快的时候是如此 ...

orca的NEVPT2梯度现在有没有人在做，以前论坛上说打算做了
作者
Author: Mikasa    时间: 2022-7-15 09:43
biogon 发表于 2022-7-15 09:39
我看22jul08的手册还说没有啊

可能只是手册没更新，我们已经用上了~
作者
Author: wzkchem5    时间: 2022-7-15 15:00
biogon 发表于 2022-7-15 02:40
orca的NEVPT2梯度现在有没有人在做，以前论坛上说打算做了

据我所知目前没有人在做。我觉得主要是NEVPT2梯度的理论比较复杂，又已经有人做过了，费半天劲做出来也不能发文章。。。DLPNO-CCSD(T)梯度虽然很可能更复杂，但起码可以发文章
作者
Author: biogon    时间: 2022-7-15 15:40
wzkchem5 发表于 2022-7-15 15:00
据我所知目前没有人在做。我觉得主要是NEVPT2梯度的理论比较复杂，又已经有人做过了，费半天劲做出来也不 ...

19年就有人发文章了，不过到现在都还没程序实现，这功能还是挺有用的
作者
Author: zjxitcc    时间: 2022-7-15 15:45
wzkchem5 发表于 2022-7-15 15:00
据我所知目前没有人在做。我觉得主要是NEVPT2梯度的理论比较复杂，又已经有人做过了，费半天劲做出来也不 ...

这是不是得怪Jae Woo Park几乎把NEVPT2的每个变种方法的解析导数都给实现and发表了。。。
作者
Author: biogon    时间: 2022-7-15 16:05
zjxitcc 发表于 2022-7-15 15:45
这是不是得怪Jae Woo Park几乎把NEVPT2的每个变种方法的解析导数都给实现and发表了。。。

那个韩国人把那些都实现了，障碍就剩编程了吧
作者
Author: zjxitcc    时间: 2022-7-15 17:01
biogon 发表于 2022-7-15 16:05
那个韩国人把那些都实现了，障碍就剩编程了吧

写对代码也很难啊
作者
Author: mfdsrax2    时间: 2022-7-15 17:03
zjxitcc 发表于 2022-7-14 11:11
2L里有说“在此之前会实现canonical CCSD(T)解析梯度”，这就是正则CCSD(T)，或者说传统CCSD(T)。传统CCS ...

我不太懂这个，请问一下像CCSD(T)这种理论，是怎么样一种函数关系式，能够写成通常意义上的f(x,y,z)吗？为什么求导如此复杂
作者
Author: wzkchem5    时间: 2022-7-15 17:03
biogon 发表于 2022-7-15 09:05
那个韩国人把那些都实现了，障碍就剩编程了吧

而且我记得他只测过很小的活性空间，类似CAS(6,5)这样的，搞不好活性空间稍微大一点就做不动了
作者
Author: zjxitcc    时间: 2022-7-15 17:08
mfdsrax2 发表于 2022-7-15 17:03
我不太懂这个，请问一下像CCSD(T)这种理论，是怎么样一种函数关系式，能够写成通常意义上的f(x,y,z)吗？ ...

你可以简单理解成100个f(x,y,z)的和 再嵌套五层复合函数。
作者
Author: beefly    时间: 2022-7-15 17:09
biogon 发表于 2022-7-15 09:39
我看22jul08的手册还说没有啊

github上有代码，但是还没并入主干版本
作者
Author: mfdsrax2    时间: 2022-7-15 17:18
zjxitcc 发表于 2022-7-15 17:08
你可以简单理解成100个f(x,y,z)的和 再嵌套五层复合函数。

如果是这样，那求导之后岂不是更复杂，如果不优化的话，算解析导数的耗时都比数值导数还要高吧，那意义在哪里，仅仅是精度高吗
作者
Author: zjxitcc    时间: 2022-7-15 17:26
mfdsrax2 发表于 2022-7-15 17:18
如果是这样，那求导之后岂不是更复杂，如果不优化的话，算解析导数的耗时都比数值导数还要高吧，那意义在 ...

解析导数形式很复杂，但计算量比数值导数小，所以只能硬着头皮做出解析导数。
作者
Author: hebrewsnabla    时间: 2022-7-15 17:30
mfdsrax2 发表于 2022-7-15 17:18
如果是这样，那求导之后岂不是更复杂，如果不优化的话，算解析导数的耗时都比数值导数还要高吧，那意义在 ...

并不是公式越复杂就越慢。一般来说，近似做得越多，公式越复杂，但是算得越快。这对于能量和解析导数都适用。
作者
Author: wzkchem5    时间: 2022-7-15 18:05
mfdsrax2 发表于 2022-7-15 10:03
我不太懂这个，请问一下像CCSD(T)这种理论，是怎么样一种函数关系式，能够写成通常意义上的f(x,y,z)吗？ ...

准确来说，不仅有多级函数嵌套，还有诸如“f(x,y)是g(x,y,z)当固定x,y只变化z时的极小值”，或者"f(x,y)是g(x,y,f(x,y))=0的解"这样的隐函数在里面，所以还涉及对隐函数求导。再者，用链式法则机械性地推出来的公式往往很慢，需要用一些技巧降低项数。
作者
Author: wzkchem5    时间: 2022-7-15 18:07
mfdsrax2 发表于 2022-7-15 10:18
如果是这样，那求导之后岂不是更复杂，如果不优化的话，算解析导数的耗时都比数值导数还要高吧，那意义在 ...

数值导数要把一个不那么复杂的函数算几十次（准确来说是原子数的6倍那么多次），解析导数只要把一个复杂的函数算一次。所以数值导数代码好写、公式好推，但（起码对于不是太小的分子来说）算得慢。
作者
Author: ionexchangeC    时间: 2022-8-10 22:39
用DLPNO-CCSD(T)-M06-2X/cc-pVTZ（M06-2X为参考态的DLPNO-CCSD(T)）按照这个方法优化，最后的结果是长碳碳键1.353埃，短碳碳键1.237埃。最后一次的ORCA输出文件和优化的过程如下。 (, 下载次数 Times of downloads: 3) (, 下载次数 Times of downloads: 3)
作者
Author: zjxitcc    时间: 2022-8-10 23:09
ionexchangeC 发表于 2022-8-10 22:39
用DLPNO-CCSD(T)-M06-2X/cc-pVTZ（M06-2X为参考态的DLPNO-CCSD(T)）按照这个方法优化，最后的结果是长碳碳 ...

为啥不用基于自旋极化UHF的DLPNO-UCCSD(T)。。
作者
Author: sobereva    时间: 2022-8-11 07:09
zjxitcc 发表于 2022-8-10 23:09
为啥不用基于自旋极化UHF的DLPNO-UCCSD(T)。。

DLPNO-UCCSD(T)比DLPNO-CCSD(T)昂贵甚巨
而且UHF错误地高估了这个体系的对称破缺特征，用U还不如用R
另外U和R我实际用CCSD优化的结果基本没差别
作者
Author: ionexchangeC    时间: 2022-8-15 11:48
又用CASSCF和NEVPT2优化了18碳环，取的活性空间是（8，8）。
CASSCF：长键1.356埃，短键1.199埃
NEVPT2：长键等于短键，1.289埃
(, 下载次数 Times of downloads: 4) (, 下载次数 Times of downloads: 7)
作者
Author: wjc404    时间: 2022-10-25 11:44
本帖最后由 wjc404 于 2022-10-25 15:34 编辑

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

另外想请教一下，现在有程序支持“DFT的基态上做CC方法”的解析梯度吗？
作者
Author: ionexchangeC    时间: 2022-11-4 23:15
wjc404 发表于 2022-10-25 11:44
是否可以用CCSD(T)-F12(RI)/cc-pVDZ-F12试一下优化(看情况启用DLPNO；这个级别的精度相当于不带F12时的TZ基 ...

按理来说，你这个计算级别就算加了DLPNO也比普通地用DLPNO-CCSD(T)/cc-pVTZ贵吧
作者
Author: wjc404    时间: 2022-11-7 23:47
本帖最后由 wjc404 于 2022-11-8 04:10 编辑
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）的活性空间可能偏小了。但再往大就比较昂贵了。
作者
Author: ionexchangeC    时间: 2022-11-8 11:38
(, 下载次数 Times of downloads: 3)
作者
Author: ionexchangeC    时间: 2022-11-8 11:39
wjc404 发表于 2022-11-7 23:47
看了一下（R1=1.22 Å，R2=1.34 Å时） RHF/def2TZVP的轨道hessian的（共3个）本征值小于零对应 ...

已上传CASSCF波函数，欢迎讨论
作者
Author: wjc404    时间: 2022-12-5 10:16
本帖最后由 wjc404 于 2022-12-5 10:24 编辑
ionexchangeC 发表于 2022-11-8 11:39
已上传CASSCF波函数，欢迎讨论
谢谢。我用ORCA跑cas(8,8)出来的波函数和这个差不多，能够收敛，但是后续的MCRPA计算发现casscf波函数不稳定（以及nevpt2计算中警告Koopmans矩阵存在负的本征值（结果未上传））。
计算结果: https://pan.baidu.com/s/1mL_uUeRGpgR_oknlpXv9rg 提取码: ari5
作者
Author: Kalinite    时间: 2022-12-19 11:01
本帖最后由 Kalinite 于 2022-12-19 11:04 编辑
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 ---
复制代码


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

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

我是为了测试超算做的这个任务，而我现在甚至忘记了测试账号的密码
作者
Author: Eclair    时间: 2024-2-10 00:14
mfdsrax2 发表于 2022-7-14 10:34
这个要写两年吗

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3