如何用Gauss直接求取并输出键长/键角的一阶二阶导数

ZetaFunction

请问各位熟悉Gauss的大佬，Gauss是否可以直接求取能量关于结构中某个键长/键角的一阶导数和二阶导数并输出呢？能量计算用的是DFT方法，导数是有解析表达式的，而且求一阶和二阶导数是Opt和Freq的基础，程序内部肯定是可以算的，我现在就是想要输出这个结果，比如我想知道内坐标下能量关于4号原子和1号原子的距离的一阶和二阶导数的值具体是多少。

我翻阅了Gauss官网的手册，上面说Opt任务的输出文件会给出每一步的一阶导数，但是我现在需要的是一个单点结构的导数，并不想做优化(否则所有一阶导数就都变成0了)。手册上说Freq会输出二阶导数，但是我检查输出文件发现是关于笛卡尔坐标的导数，虽然也可以手动转换但还是有点麻烦最好可以直接输出。手册上还说在输入文件末尾写上“B 4 1 D”就会计算4号原子和5号原子距离的二阶导数，但是我照做后翻阅了整个输出文件都未见计算结果。

1. %nprocshared=40
   
2. %mem=160GB
   
3. %chk=BD-Der.chk
   
4. # M062X/gen geom=modredundant
   
5. 
   
6. Der
   
7. 
   
8. -1 1
   
9. C                  0.71241700   -0.23209900    0.35036600
   
10. O                  1.21233200   -1.20834600   -0.22999200
    
11. H                  0.54006600   -0.28116200    1.44293400
    
12. C                 -1.27239600   -0.05182500    0.04828100
    
13. N                 -2.42123100    0.00930900   -0.09903900
    
14. C                  1.08316200    1.17684700   -0.10635200
    
15. H                  0.48030600    1.93470600    0.39556400
    
16. H                  2.14204800    1.33286000    0.12512100
    
17. H                  0.94843600    1.25766400   -1.18418300
    
18. 
    
19. B 4 1 D
    
20. A 4 1 2 D
    
21. D 3 1 2 4 D
    
22. A 5 4 1 D
    
23. A 4 1 6 D
    
24. 
    
25. @ma-TZVP.txt

复制代码

希望有大佬能帮我解答，如果Gauss确实做不到，有没有其他的量化程序可以实现这样的功能？

mizu-bai

force nosymm 可以满足你的需求

sobereva

不要把Gaussian写成Gauss

用freq和#P可输出，参考北京科音基础量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）的ppt

ZetaFunction

sobereva 发表于 2023-8-19 05:01
不要把Gaussian写成Gauss

用freq和#P可输出，参考北京科音基础量子化学培训班（http://www.keinsci.com/ ...

感谢Sobereva老师的解答，我在输入文件写上“#P freq”后逐行检查输出的log文件，还是只能找到笛卡尔坐标下的海森矩阵和一阶导数，并未找到内坐标下定义的海森矩阵。另外海森矩阵的单位是Hartree/Bohr^2，请问如何转换到Hartree/A^2？

ahxb

输入文件中写上freq=modredundant后，chk文件保存了内坐标下的受力和海森矩阵。计算完毕后使用formchk将chk文件转化为fchk文件，可在fchk文件中看到Internal Forces和Internal Force Constants字段（以下三角矩阵形式表示，可与输出文件中的结果比对，但有效数字更多）。冗余内坐标定义可在输出文件开头附近看到，fchk文件中也有Redundant internal coordinate indices和Redundant internal coordinate字段，前者是冗余内坐标定义，四个一组；后者是冗余内坐标数值。
也可以执行rwfdump ***.chk 511.txt 511R，在文本文件511.txt末尾会给出一个数组，第一个数是SCF收敛后的能量；第四个数开始是在内坐标中的受力，长度为冗余内坐标个数NVAR；然后是内坐标下的二阶力常数（以下三角矩阵形式表示），长度为NVAR*(NVAR+1)/2。你要是还想获得更高精度的结果，还可以从rwfdump输出的文本文件中Number为511的部分读取Base和End，然后自行写程序读取二进制文件中从Base字节到End字节的部分，与文件最后给出的结果完全相同，是Gaussian计算出的原始结果，精度只受双精度浮点数本身的限制，一般可达15~16位有效数字。
1 Hartree/Bohr^2 = 1 Hartree/(0.529 angstrom)^2 = 3.57 Hartree/angstrom^2，想要更精确的转换关系可以去查Bohr的高精度定义。

sobereva

ZetaFunction 发表于 2023-8-19 19:59
感谢Sobereva老师的解答，我在输入文件写上“#P freq”后逐行检查输出的log文件，还是只能找到笛卡尔坐标 ...

拿具体输入输出文件说事。ppt里该说的都说了