如何计算范德瓦尔斯力？ 第1页

本文已经升级，详细讨论请移步专栏。

首先，这个问题不用量子力学也能解释。（为讨论问题方便，以下全用Gauss单位制，部分地方用了原子单位制）

根据电动力学的多极展开理论，两个电中性 分子（比如惰性气体这种单原子分子）可以分别展开成多极矩，领头的非零项就是四极－四极相互作用。

（这里需要注意模型，必须使用电荷连续分布的电中性带电体模型，不能把原子核和电子单独考虑。否则必然出现偶极－偶极相互作用，导致高估相互作用能。）

咱们先不用多极展开的一般公式，暴力硬算一个简单的情况：

不考虑它们之间的角度，用准一维的模型如上图，八个单位电荷如上图，可以直观的写出两个四极矩的相互作用能：

左边左上角正电荷跟右边四个电荷相互作用：

左边右上角负电荷跟右边四个电荷相互作用：

根据对称性下边两个电荷跟右边四个电荷的相互作用之和等于上边八项和。

总的相互作用能：

将所有的分式用泰勒级数展开，就可以得到：

因为此处设

代入 ，别忘了“1”是 ，故有（如果直接设两分子间距为 更不容易出错）：

两个四极矩的相互作用力就有：

这就是从静电学看出来的分子间相互作用力。

如果用严格的多极展开，可以得到（A. Towns, The theory of intermolecular forces, Oxford University Press, 2013）：

上式用了SI单位制，分子上两项是A、B分子的标量四极矩。这个表达式也可以参考Jackson的书第四章内容，特别是公式（4.10）、（4.21）和（4.24），自己算出来。

量子力学的讨论更简单一些，只需用到偶极－偶极相互作用。此时最简单模型中共有四个电荷，分别属于两个原子核和两个电子。应用Born-Oppenheimer近似的思想，原子核质量为电子的1823倍，近似不动。体系中是两个量子的偶极子相互作用（以下参考Cohen-Tannoudji，附录 ）。

偶极－偶极相互作用能由如下方法得到：

偶极子A在距离 处产生静电势：

电场可求：

这里 为 的方向矢量。

上式的计算用到了：

使用定态微扰法，体系哈密顿为：

这里对上述偶极偶极相互作用势可以假设 沿着 轴，进一步化简为：

（复杂的情况可以引入三个角度，然后按第二段的方法用微扰展开求统计平均）

零级哈密顿的本征方程为：

考虑两个 态的氢原子，则一阶修正项为 。这是因为各坐标分量的算符平均值为 。

二阶修正为：

这里" "代表求和中不包含基态。而且可以看出，因为 ，分母恒负，而分子又是一个模方，恒大于等于零，这个二阶修正项取负值。即：

经典电动力学和量子力学微扰法的结果正好差一个 。但是需要注意，上述讨论仅限于两个分子距离较远，但是又不能太远的情况。离得过近了显然会出现排斥（不远不近可能会有轨道相互作用，比如本人和同事研究过的卤原子－水三电子化学键 JCP 138 (14), 141102），离得远了则有所谓延迟效应存在，就是A分子在B分子处产生的场被B分子内电子的运动扰动，这个扰动要以有限的速度即光速传回A，改变A在B处产生的场。根据细致的讨论，会发现这个距离一般不能超过100 nm（详见Cohen-Tannoudji）。

在实际的应用中，这个吸引势一般取为 的形式，如Lonnard-Jones势，但是两个中性分子靠近后会排斥，需要添加排斥势，这个势无法通过普通的量子力学计算得到。一般猜测具有 的形式。近年也有用三项和（Kenneth Merz组2014年的工作）：

的形式来取得更精确的近似。

详细的讨论需要上量子场论，请参考近年来Tkatchenko的工作。Alexandre TKATCHENKO

结束语

实际上，经典物理可以帮我们建立许多体系的基本模型，虽然会存在各种误差，但是正如Weisskopf说的：“物理学就是列车时刻表。列车总会误点的，但是只有你拿着列车时刻表，才知道它误点了啊！”拿到问题先使用经典理论分析一番，一般就能获得一个问题的概貌了。

再比如说，氢原子问题可以用匀速圆周运动＋库伦吸引力＋角动量量子化分析一波，可以获得波尔半径、基态能量等重要信息，跟量子力学的结果相差无几。再比如氢分子离子，也可以用库伦力分析得到其基本的成键信息。近年量子光学大佬M. O. Scully还使用Bohr-Sommerfeld量子化计算得到了一系列正确的基态双原子分子。而计算电子结构时，先上DFTB、LDA之类估计一下，估算大范围内的结果，再在特别感兴趣的地方加细，也是我们常用的技巧。

复杂的理论往往是用来定量，但是算多了总会迷失在数字中，忘记物理本质。个人因为数学不是特别好，总喜欢用微扰展开硬算一些简单的情况（1维、2个粒子等等）。但经验告诉我，这样往往可以把握住问题的大方向。老乡（虽说他不认，但是他的确是在长春出生并上学到小学毕业）彭桓武也总说，1，2，3为无穷大；遇到问题先估算等等。

高等量子相声学老师也曾说，很多公式fancy的文章，往往其作者会熬夜算微扰展开到20阶，在凌晨五点得到最终结果后，换上一套路径积分之类的高级公式投稿，把微扰计算的草稿纸严严实实地锁在抽屉里。

范德瓦尔斯力，又被称作分子间作用力，是一种中性分子或者原子之间的相互作用力，其产生的根源依然是分子或者原子之间的静电相互作用，只不过此时发挥作用的是分子或者原子中的偶极矩。

范德瓦尔斯相互作用与化学键相比非常微弱，但在一些分子晶体当中，由于分子与分子之间存在范德瓦尔斯相互作用，因此需要考虑其对结构所造成的影响。此外，在对层状结构或者用两个单层材料制作而成的异质结做第一性原理计算时一般也要考虑范德瓦尔斯相互作用。对于单层的体系，一般可以不用考虑范德瓦尔斯。

下面我们就将介绍VASP中几种常见的加范德瓦尔斯相互作用的方法；

1. DFT-D2
   该方法用法简单，只需在INCAR中加入开关IVDW=1（或者IVDW=10）即可。但是这里需要注意两点：(i)该方法只适用于5.2.11以后的版本；(ii)在一些比较老的版本中，VASP只给出了从氢（元素符号：H）到氙（元素符号：Xe）的范德瓦尔斯修正相关参数，本人亲测5.3.3版本的VASP在加上IVDW=1的开关计算重元素时，VASP会提示需要用VDW_C6来为C6取值，并且计算终止，但采用5.4.1版本的VASP不会出现该提示，软件会正常运行。
   
2. DFT-D3
   该方法只需在INCAR中加入开关IVDW=11即可。该方法是DFT-D2的增强版，5.3.4及以后的版本可用。
   

在这里大家注意，我们并没有在INCAR里用到GGA开关，这表明我们并没有改变交换关联势，以下四种范德瓦尔斯修正就和上面的两个不同，需要用到GGA开关。

3. vdW-DF

采用该范德瓦尔斯修正需在INCAR中加上以下开关：

GGA=RE

LUSE_VDW=.TRUE.

AGGAC=0.0000

4. optB88-vdW

采用该范德瓦尔斯修正需在INCAR中加上以下开关：

GGA=BO

PARAM1=0.1833333333

PARAM2=0.2200000000

LUSE_VDW=.TRUE.

AGGAC=0.0000

5. optB86b-vdW

采用该范德瓦尔斯修正需在INCAR中加上以下开关：

GGA=MK

PARAM1=0.1234

PARAM2=1.0000

LUSE_VDW=.TRUE.

AGGAC=0.0000

6. vdW-DF2

采用该范德瓦尔斯修正需在INCAR中加上以下开关：

GGA=ML

LUSE_VDW=.TRUE.

Zab_vdW=-1.8867

AGGAC=0.0000

注意，vdW-DF、optB88-vdW 、optB86b-vdW和vdW-DF2在使用过程中需要在所在的文件夹下放置一个名叫vdw_kernel.bindat的文件，该文件我已经上传到QQ群“凝聚态与VASP使用指南”的群文件中，群号为427494979，需要的同学可自行下载。

此外，本人测试的结果表明，在考虑DFT-D2或者DFT-D3计算的过程中可同时考虑自旋轨道耦合（SOC），但是在vdW-DF、optB88-vdW 、optB86b-vdW和vdW-DF2中无法同时考虑自旋轨道耦合进行计算（VASP倒是运行了，但是并没有按照预期出现修正）。有没有同学遇到过范德瓦尔斯+SOC的问题？对这方面比较明白的同学可以将答案发在留言板，欢迎大家交流。