求助DFT-D3能量、解析导数和解析Hessian如何写表达式

ldatea

最早写的（分割线后面的段落）关于G09 D.01 DFT-D3校正能的叙述是错的，G09 D.01 DFT-D3的能量和梯度计算是没有问题的。被我搞出问题来是因为我为了算DSD-PBEP986-D3(BJ)把环境变量永久改了。出现这么蠢的问题是因为我想当然了。
按照sob老师的做法修改环境变量。才能够正常算DSD-PBEP986-D3(BJ)，但是修改环境变量以后，使用g09w.exe运行还是会出错，出错原因是D3参数没有定义
我就想会不会是g09w.exe会自动把环境变量重新设一遍，之前这个问题让我很困扰。
不过使用g09以后就不存在这个问题了。
最终发现这个问题还是因为，我最早算的一批数据里使用了DSD-PBEP986-D3(BJ)，我原本以为g09 D.01的D3能量也是错的，于是重新用g16算了一遍。结果，两者的结果一致。
到这个时候我就明白了，肯定是我之前设的永久生效的环境变量出事了。
    （以下是自己尝试解决问题的过程，与求助内容关系不大）
——————————————————————————————————————————————————————————————————     之前看了sob老师的《DFT-D色散校正的使用》http://sobereva.com/210，发现自己只有G09D.01版是可以正常使用的，里面给的一个external的示例，在几何优化时使用XDM校正，但是网页已经打不开了。     今天尝试给G09 D.01使用外部程序提供的DFT-D3Hessian，结果发现没有现成的程序，Grimme的DFT-D3程序只能获得能量和一阶解析导数。然后尝试借助ORCA的DFT-D3给Gaussian传数据，试了一下sob老师的gau_orca接口，发现DFT-D3是内嵌的，算导数的时候是scf导数和DFT-D3导数合在一个程序里算的，这样就没办法单独获得D3部分的导数，如果非要这么干，就只能全部拿ORCA计算。
      接着就开始研究Grimme的DFT-D3程序，将DFT-D3程序输出的部分数据用Excel手动计算一下，看看这个公式到底是怎么使用的。
      DFT-D3的公式其实不难理解，但是我没看到里面用的单位，就打算用Excel来确定里面数据的单位（特别是RAB）。
最后的结果很明显，原子间距用的是a.u.不是gaussian默认的埃。其实这个点很明显，因为Grimme的DFT-D3程序虽然输入坐标的时候单位是埃，但会转化为a.u.再计算，这在输出里就有体现。然后是能量公式里的1/2，我发现按原文带1/2的公式算，和最终的输出相比，能量少了一半。附件里面是两个最简单模型的计算，分别是相距1埃和2埃的两个氢原子。DFT-D3的输出都放在“DFT-D3程序的输出”文件夹下，DFT-D3自动生成的输出文件为了区分已经添加了HH1，HH2以区分。
所有计算的泛函都是B3LYP，零阻尼的结果都带zero,没标明的都是BJ阻尼。


解析导数我瞎算的，因为把C6AB和C8AB当常数了，但比较HH1.txt和HH2.txt两个文件可以很明显的看出C6AB和C8AB是和原子间距有关的。
然后分别用g16C01和g09D.01算了一下，发现G16不管加不加G09default，计算的数据和DFT-D3完全一致。但是G09的数据不管是BJ阻尼还是零阻尼全是对不上的。我把G16的数据减掉G09的数据，也没看出来G09的计算结果是怎么来的。
目前在Linux系统上使用external很顺利，不过很迷的是某瓜自带的脚本extuffex输出的格式居然不是标准格式（对于能量和偶极矩，手册上的标准格式是4D20.12）

1. awk '{if($1=="Energy=") printf "%20.12e%20.12e%20.12e%20.12e\n", $2,0.0,0.0,0.0}' tempuff2 > tempuff3

复制代码

其实也是大同小异，只不过把D换成了E，而且个位可以不为0，不过awk命令只“认识”“%20.12e”格式，如果你使用nD20.12，转换后指数就是0了。
不过DFT-D3的xyz输入文件也可以使用带E的这种格式，sob老师的gau_xtb也用的是这种。看起来这两种格式似乎是通用的。

但是Windows系统非常坑，D.01自带的extuffex.bat直接运行会报错，而且报错很怪，是L402不能打开 output file，但实际上是 output file只含能量。（我不清楚是不是我操作有误）
不过稍微改一下代码，这一步就没问题了，但是接下来L716又会报错，我很懵逼。

gaussian还可以用iop(7/32=x)输出能量，坐标，力，力常数等一系列值。x=2会在L716报错，而且这个报错很奇怪，输出文件里什么也没有，Linux下终端会出现

1. FIO-F-246/formatted write/unit=7/infinite format scan for edit descriptor.
   
2. File name = 'fort.7',    formatted, sequential access   record = 4
   
3. In source file frcout.f, at line number 1604
   
4.   /home/wyn/apps/g16/l716.exe() [0x188f517]
   
5.   /home/wyn/apps/g16/l716.exe() [0x18492bf]
   
6.   /home/wyn/apps/g16/l716.exe() [0x1826b2f]
   
7.   /home/wyn/apps/g16/l716.exe() [0x182ddc6]
   
8.   /home/wyn/apps/g16/l716.exe() [0x182db1e]
   
9.   /home/wyn/apps/g16/l716.exe() [0x187c4fe]
   
10.   /home/wyn/apps/g16/l716.exe() [0x188db20]
    
11.   /home/wyn/apps/g16/l716.exe() [0x187c662]
    
12.   /home/wyn/apps/g16/l716.exe() [0x1865265]
    
13.   /home/wyn/apps/g16/l716.exe() [0x41b3d9]
    
14.   /home/wyn/apps/g16/l716.exe() [0x411fa2]
    
15.   /home/wyn/apps/g16/l716.exe() [0x411ea4]
    
16.   /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf3) [0x7f399af270b3]
    
17.   /home/wyn/apps/g16/l716.exe() [0x40c5e9]

复制代码

看起来像是FileIO程序出了问题。
不过用x=3就可以了，L716会产生一个PES.out文件。

Daniel_Arndt

根据我个人记忆，Gaussian 09 D.01里面的D3可能是有点问题的，但是Gaussian 09 E.01据说修复了这个bug。楼主估计是用的Gaussian 09 D.01，没有用E.01。