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

ionexchangeC · 发表于 Post on 2022-7-13 23:29:21

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

18C.f90 (51.76 KB, 下载次数 Times of downloads: 9)

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

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

复制代码

最后程序的输出是

-0.685090685940D+03 0.134360000000D+01 0.122420000000D+01
-0.685088824064D+03 0.135360000000D+01 0.122420000000D+01
-0.685091935687D+03 0.134360000000D+01 0.123420000000D+01
-0.685091376254D+03 0.133360000000D+01 0.122420000000D+01
-0.685086911391D+03 0.134360000000D+01 0.121420000000D+01
-0.685093443649D+03 0.133360049157D+01 0.123410084787D+01
-0.685094798943D+03 0.132032800904D+01 0.124176231458D+01
-0.685094995375D+03 0.131259737079D+01 0.124433184094D+01
-0.685095034243D+03 0.131324281292D+01 0.125431098942D+01
-0.685095381505D+03 0.130901412440D+01 0.125156104932D+01
-0.685095043825D+03 0.130676166617D+01 0.126130406897D+01
-0.685095617327D+03 0.130113934072D+01 0.125772447223D+01
-0.685096084259D+03 0.128547177783D+01 0.126955992432D+01
-0.685094926437D+03 0.125415034613D+01 0.129324894559D+01
-0.685094328761D+03 0.126621797661D+01 0.127341240296D+01
-0.685095912811D+03 0.128405212975D+01 0.127945864137D+01
-0.685096232455D+03 0.127922217111D+01 0.127736648668D+01
-0.685096136254D+03 0.127312956531D+01 0.128529618747D+01
New RHO = 1.0000D-03 Current number of function evaluations = 18
Least value of F = -6.850962324546280D+02
The corresponding X array is:
1.279222D+00 1.277366D+00
-0.685096233410D+03 0.127758652046D+01 0.127882988316D+01
-0.685096234289D+03 0.127809725865D+01 0.127837670778D+01
-0.685096220707D+03 0.127736049474D+01 0.127770055769D+01
-0.685096223664D+03 0.127897725857D+01 0.127885168162D+01
New RHO = 1.0000D-04 Current number of function evaluations = 22
Least value of F = -6.850962342889070D+02
The corresponding X array is:
1.278097D+00 1.278377D+00
-0.685096234359D+03 0.127830643932D+01 0.127823790676D+01
-0.685096234358D+03 0.127837923057D+01 0.127816933975D+01
-0.685096234369D+03 0.127823904017D+01 0.127816403285D+01
-0.685096233929D+03 0.127816395443D+01 0.127809798641D+01
-0.685096234265D+03 0.127830070069D+01 0.127817211338D+01
-0.685096234304D+03 0.127816288187D+01 0.127822883960D+01
-0.685096234345D+03 0.127830944408D+01 0.127816255353D+01
At the return from NEWUOA Total times of function evaluations = 29
Least value of F = -6.850962343688360D+02
The corresponding X array is:
1.278239D+00 1.278164D+00

复制代码

wzkchem5 · 发表于 Post on 2022-7-14 01:34:54

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

sobereva · 发表于 Post on 2022-7-14 06:29:37

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

hebrewsnabla · 发表于 Post on 2022-7-14 09:56:12

本帖最后由 hebrewsnabla 于 2022-7-14 10:06 编辑

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

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

mfdsrax2 · 发表于 Post on 2022-7-14 10:34:52

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

这个要写两年吗

zjxitcc · 发表于 Post on 2022-7-14 10:52:24

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

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

mfdsrax2 · 发表于 Post on 2022-7-14 10:54:16

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

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

zjxitcc · 发表于 Post on 2022-7-14 11:00:21

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

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

mfdsrax2 · 发表于 Post on 2022-7-14 11:07:01

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

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

zjxitcc · 发表于 Post on 2022-7-14 11:11:02

本帖最后由 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)解析导数。

wzkchem5 · 发表于 Post on 2022-7-14 17:00:16

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

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

Mikasa · 发表于 Post on 2022-7-15 09:10:26

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

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

biogon · 发表于 Post on 2022-7-15 09:39:24

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

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

biogon · 发表于 Post on 2022-7-15 09:40:51

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

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

Mikasa · 发表于 Post on 2022-7-15 09:43:37

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

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

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[算法与编程] 一个用DLPNO-CCSD(T)优化18碳环的结构的思路

评分 Rate

评分 Rate