分子动力学模拟是化学吗？ 第1页

分子动力学忽略电子，而我们常说的化学反应的本质是电子的变化，那分子动力学的研究和化学有什么关系呢？

题主的全部疑惑可能来源于对分子动力学（MD）的错误理解：MD 并没有忽略电子！

我可以举一个最简单的化学反应的例子：氢分子解离成两个氢原子（ ）。

我们考虑这个过程中各种能量的变化。

首先是两个氢原子核的相互作用。这是完全符合题主要求的“忽略电子”。会得到如下图所示的曲线：随键长增大能量逐渐下降，这是因为没有电子的时候两个带正电荷的核之间只有库伦斥力。。。距离越大斥力越小。。。这时候完全没有 chemistry（没有 chemical bond！）。

接下来我们加入电子，考虑“两个氢原子核 + 两个电子”这一体系的总能量。对比上图我们可以看到总能量随键长的变化曲线上在 的时候出现了一个极小值。这说明在这个键长时氢分子是稳定的。你问我这是什么。。。这当然就是所谓的“化学键“。。。所以正如题主所说的，加入电子之后就有了“化学”。

上面第二张图虽然有化学键，但是那是因为我在背后帮你解了“薛定谔方程”（代码附在后面了）。现在回到分子动力学（假设是经典 MD，量子情况可以看 @蘭一 的回答），你要问的问题是，假设我不想/不会/懒得解薛定谔方程，那么我该如何描述好氢分子的解离（也就是你说的“化学”）？

答案是我可以拟合啊！！！这不就是个单变量函数嘛，我都学过 machine learning 了这还难得倒我。。。（ 狗头）

随手 google 一下你会搜到一个叫 Morse potential 的东西：

它是前人用来拟合这种解离曲线的一个比较合理的函数形式：

(Eqn 1)

作为一个出色的数据科学家（雾），我迅速修改了我的代码使它能够拟合 Morse potential，然后我得到了下图的结果：

惊喜不惊喜？意外不意外？Morse potential 非常准确地描述了氢分子解离过程中能量的变化。换句话说，它准确地描述了化学键。

现在你仔细想一下，如果我给你了我拟合出来的 Morse potential 的参数，你只需要氢分子的键长 就能准确算出它的总能量从而能描述和它有关的“化学过程”。最重要的是，这个过程中你并不需要解薛定谔方程！所以看上去，你忽略了电子。但其实电子的信息早已经在 Morse potential 的参数里面了。

这里拟合的 Morse potential 可以看成描述氢分子的最最最最最最最简单的分子动力学“力场”。有了它，你就可以在波恩-奥本海默近似下描述氢分子中两个原子的运动了。

对更复杂的体系，显然没有像上面描述的这么简单的函数形式来拟合多体薛定谔方程的解，不过不管怎么说人们还是搞出了很多有用的“力场”使得基于它们进行的分子动力学模拟能够给出有意义的化学。

所以总的来说，MD 是可以用来描述化学的。

====================================

       """ 'Exact' quantum mechanical treatment of H2 dissociation in minimal basis. The results are fitted using a Morse potential. """  import numpy as np import scipy.optimize as scopt  from pyscf import gto, scf, fci  if __name__ == "__main__":     basis = "sto-3g"     bL_list = np.arange(0.4,3.01,0.1)     e_tots = np.zeros(bL_list.size)     e_nucs = np.zeros(bL_list.size)      # loop over various bond lengths     for i,bL in enumerate(bL_list):         geom = "H 0 0 0; H 0 0 %f" % bL         mol = gto.M(atom=geom, basis=basis, spin=0)         mf = scf.RHF(mol)         mf.kernel()         mc = fci.FCI(mf)         mc.kernel()          e_tots[i] = mc.e_tot         e_nucs[i] = mf.energy_nuc()      # fit a Morse potential!     def morse_pot(rs, De, alp, re, E0):         return De * (1. - np.exp(-alp * (rs - re)))**2. + E0      p0 = [0.2, 1.0, 0.8, -1.1]     res = scopt.curve_fit(morse_pot, bL_list, e_tots, p0=p0)     e_tots_morse = morse_pot(bL_list, *res[0])      # plot! (略。。。)