非简谐振动下热矫正量的计算方式?

wuzhiyi

最近刚入坑量子化学，想请教一下关于氢键系统自由能的问题，对于一个普通的氢键系统比如D-H···A。我们计算的这个系统自由能的方式是将结构优化后进行频率计算算热矫正量，然后将优化的结构算单点，体系的自由能=热矫正量+单点能。
我对于这种算法的理解是，我们可以画出D-H···A中H在D和A之间的势能平面。X轴是H在D>A之间的位置(x)，Y轴是系统的势能(H)。当H靠D非常近的时候，系统的势能很大，等H离开D接近1A的时候，是系统的势能最小值，然后随着H逐渐远离D，因为D-H键的断裂，势能开始增加。然后随着H逐渐靠近A，随着H-A键的形成，系统势能(H)下降，但因为系统是D-H···A，所以在靠近A处的势能最小点比靠近D处的势能最小点要高上不少。
在计算体系D-H···A自由能的时候，我们定义G=∫exp(-H(x)/kbT)dx，同时将系统H(x)简化为H(x)=C+k(x-x0)^2。期中C就是黄色箭头势能最低点处的势能，相当于优化之后系统的单点能。而k(x-x0)^2则是一个harmonic oscillator，相当于将蓝色箭头处的振动纳入考虑，其中x0是黄色箭头势能最低点处的x值，而spring strength则可以通过二阶导数得到频率计算，其体现就是热矫正量。

对于普通的氢键系统，这样的计算问题不大，如下图所示，因为如果画出势能平面的话，确实一个harmonic oscillator可以比较准确的描述系统的势能平面。虽然系统有两个局部最小点，但因为两个势能最小点之间有一个局部最大点，所以可以基于局部最大点将系统分为D-H···A和D···H-A两种情况分别讨论。


但我最近在研究的一个氢键系统比较诡异，类似low-barrier hydrogen bond，其势能平面如下图所示。一，在D-H···A和D···H-A之间不存在局部最大点，所以没法分两种情况讨论。二，系统势能平面与harmonic oscillator相去甚远，非谐振效应非常显著。

请问一下，这种情况如何计算算热矫正量。谢谢

Accelerator

最简单的方法是用Gaussian自带的VPT2方式考虑非谐性 freq=anharmonic

sobereva

如果第二个图里势阱深度很浅的话，实际被研究的温度下也会有大概率运动到右边去，即便用非谐振恐怕结果也不理想。这时候需要根据扫出来的势能曲线自己求解振动态能级，并且自行按照定义算振动配分函数，并得到振动对热力学量的贡献

wuzhiyi

sobereva 发表于 2019-6-14 13:17
如果第二个图里势阱深度很浅的话，实际被研究的温度下也会有大概率运动到右边去，即便用非谐振恐怕结果也不 ...

第二个图的势阱深度只有2kcal/mol，这几天试了试freq=anharmonic，好像要算很久的样子。
扫完了势能曲线后，这两天读了读书，发现好像要在每一个能级numerical解薛定谔方程算激发能，想问一下sob老师可不可以推荐对已知非简谐势能曲线求热校正量相关的文献？

YanLiu

最近也在看一些非谐振的文章，这两篇可能对你有帮助（虽然我还没有看明白）。顺便问一下，这个问题楼主解决了吗？

wuzhiyi

YanLiu 发表于 2019-11-16 09:59
最近也在看一些非谐振的文章，这两篇可能对你有帮助（虽然我还没有看明白）。顺便问一下，这个问题楼主解决 ...

我现在用的方法在原理上不太严格，不过我是要比大小，而不是具体量化分析，所以我们老板还算满意，就看审稿人怎么说了。
我就假设在x坐标轴上和x轴无关的震动所带来的热矫正量处处相等。这明显是错的，因为x轴正交方向上的震动和x轴的位置有关。在x轴方向上的震动用电子能量来代替，就是先减去在x轴上的简协震动的entropy，再加上根据势能平面算出来的在x轴上运动的entropy。

YanLiu

大致明白楼主的做法，谢谢楼主。